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CHAPTER  1 


1.1  Introduction 

Development  of  intelligent  material  systems  and  structures  has  been  followed  by  various 
government  agencies  and  private  corporations  in  the  last  decade.  Intelligent  structures  are 
defined  as  the  ones  which  incorporate  actuators  and  sensors  that  are  highly  integrated  into 
the  structure  and  have  structural  functionality,  as  well  as  highly  integrated  control  logic, 
signal  conditioning,  and  power  amplification  electronics.  Such  actuating,  sensing,  and 
signal  processing  elements  are  incorporated  into  a  structure  for  the  purpose  of  influencing 
its  mechanical,  thermal,  optical,  chemical,  electrical  or  magnetic  characteristics. 

Several  technical  developments  have  been  combined  to  establish  the  potential  feasibiUty 
of  intelligent  structures.  The  first  is  a  transition  to  laminate  materials.  In  the  pzist, 
structures  were  built  from  large  pieces  of  monolithic  materials  which  were  machined,  forged, 
or  formed  to  a  finsJ  structural  shape,  malcing  it  diflBcult  to  incorporate  any  active  element 
inside  the  structure.  However,  in  the  pzist  thirty  yeairs  a  transition  to  laminated  material 
technology  has  occurred.  Laminated  materials  allow  for  the  esisy  incorporation  of  active 
elements  within  the  build  up  of  structural  forms. 

Exploitation  of  the  off  diagonal  terms  in  the  material  constitutive  relations  is  a  second 
trend  which  madces  the  development  of  intelligent  structures  promising.  The  full  consti¬ 
tutive  relations  of  a  material  which  include  characterizations  of  its  mechanical,  optical, 
electromagnetic,  chemical,  physical,  and  thermeil  properties.  Normally,  resezirchers  have 
focused  only  on  block  diagonal  terms.  For  example,  those  interested  in  exploiting  a  ma¬ 
terial  for  its  structural  benefits  have  focused  only  on  the  mechanical  characterization,  and 
those  interested  in  exploiting  its  electrical  properties  have  focused  on  the  electrical  charac¬ 
terization.  However,  much  gain  can  be  made  by  exploiting  the  off  block  diagonal  terms  in 


the  constitutive  relations,  which  for  example,  couple  the  meclianical  and  electrical  proper¬ 
ties  or  couple  mechanical  and  thermal  properties.  The  characterization  and  exploitation 
of  these  off  diagonal  material  relations  had  led  to  much  of  the  progress  in  the  creation  of 
intelligent  structures. 

The  third  and  perhaps  most  obvious  advance  has  been  in  the  electrical  engineering 
and  computer  science  disciplines.  These  include  the  development  of  microelectronics,  bus 
architectures,  switching  circuitry,  and  fiber  optic  technology.  Also  central  to  the  emergence 
of  intelligent  structures  is  the  development  of  information  processing,  artificial  intelligence, 
and  control  disciplines.  The  sum  of  these  three  developing  technologies  (the  transition  to 
laminated  materials,  the  exploitation  of  the  oif-diagonal  terms  in  materieJ  constitutive 
relations,  and  the  advances  in  microelectronics)  hais  created  the  enabling  infrastructure 
based  on  which  intelligent  structures  can  develop. 

A  wide  variety  of  applications  exist  for  intelligent  structure  technologies.  These  include 
aeroelastic,  control  emd  maneuver  enhancement,  reduction  of  vibrations  and  structure 
borne  noise,  jitter  reduction  in  precision  pointing  systems,  shape  control  of  plates,  trusses 
smd  lifting  surfaces,  isolation  of  offending  machinery  euid  sensitive  instruments,  and  robotic 
control. 

For  a  review  of  the  status  of  smart  structures  and  technology  overview,  the  following  ref¬ 
erences  may  be  consulted:  Agnes  &  Silva  [1992],  ’’Aircraft  Smart  Structures  Research  in  the 
USAF  Wright  Laboratory”;  Crawly  [1992],  ’’Intelligent  Structures,  a  Technology  Overview 
and  Assessment”;  Lazarus  &  Napolitano  [1993],  ’’Smart  Structures,  an  Overview.” 

Three  main  development  eirezis  that  require  further  research  are:  (1)  basic  materials, 
(2)  component  devices,  and  (3)  modeling  techniques.  There  is  a  need  to  develop  basic 
smart  materials  with  improved  properties  and  characteristics.  Some  of  the  more  impor¬ 
tant  features  to  be  improved  include  fatigue  life,  fracture  toughness,  structural  weight, 
required  voltage,  linearity,  eind  anisotropy.  These  materials  should  be  developed  in  order 
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to  create  improved  intelligent  structure  component  devices  such  as  actuators  and  sensors. 
An  emphasis  should  be  placed  on  the  use  of  actuator  and  sensor  components  for  structur£Ll 
control.  Accurate  models  need  to  be  developed  so  that  smart  material  components  may 
be  utilized  elfectively.  Work  is  needed  in  developing  constitutive  relations  as  well  as  fiill 
non-lineax  models  to  meet  these  needs.  This  research  is  directed  toward  addressing  some 
of  these  needs  els  discussed  in  the  following  sections. 

In  an  intelligent  structure,  the  presence  of  active  elements  (actuators,  sensors,  and 
processors)  impact  the  host  structure  by  increasing  the  mass  of  the  system  and  interfer¬ 
ing  with  the  load  path  and  potentially  introducing  new  structural  discontinuities  which 
must  be  accommodated.  This  may  potentially  change  the  fatigue  and  fracture  toughness 
characteristics  of  the  host  material. 

The  complex  microscale  interactions  between  sensors,  actuators,  and  the  host  materi¬ 
als,  and  in  particular  the  inherent  material  and  geometric  nonlinearities  within  intelligent 
material  systems,  must  be  characterized  before  the  technology  can  reach  its  full  potential. 

Initial  experimentEil  investigations  have  indicated  high  strain  concentration  at  interface 
between  sensors  and  host  materials.  Czaumek  et  al.  [1988]  at  Virginia  Tech  have  used 
interference  imaging  for  a  cross-ply  graphite/epoxy  laminate,  fabricated  with  a  jacketed 
gletss  fiber  waveguide  embedded  between  and  perpendiculeu’  to  the  center  two  plies.  Strains 
of  approximately  0.05  were  indicated  at  the  fiber-matrix  boundary  for  an  applied  load  equEil 
to  half  the  failure  load  of  specimens  tested.  They  concluded  that  these  large  interface  strain 
concentrations  may  pose  significant  limits  on  the  long  term  structural  integrity  of  materials 
containing  embedded  sensor  fibers. 

Singh  H.  et  al.  [1991]  combined  Moire  interferometry  Eind  the  Fourier  transform  fringe 
interpretation  method  to  produce  high  resolution  displacement  measurement  system  for 
laminated  composites  with  embedded  sensors.  Their  findings  show  stress  concentrations 
factors  ELS  much  as  8.8  in  the  vicinity  of  fiber-host  material  interfeice.  The  largest  character- 
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istic  length  of  the  influence  was  found  to  be  about  11.5  times  fiber  diameter.  These  strain 
concentrations  at  sensor-host  material  interface  warrant  concern  regarding  micro-cracking 
due  to  various  in-use  loading  conditions. 

Other  studies,  Jenson  [1992],  have  measured  reduction  in  compression  strength  due  to 
embedment,  for  the  worst  orientation,  of  up  to  70  %  .  This  reduction  would  seem  to  be 
due  to  the  bending  of  load-bearing  longitudinal  fibers  around  a  transverse  sensory  fiber, 
creating  a  prebuckled  site. 

These  experimental  results  indicate  the  need  for  research  work  to  address  the  mechainics 
that  govern  the  complex  interactions  between  sensors,  actuators,  and  the  host  materials 
for  smart  structures.  The  research  should  consider  the  inherent  material  and  geometric 
nonlinearities  at  micro  and  macro-scales.  Fundamental  understanding  of  these  interactions 
will  provide  opportunities  for  further  developments  of  intelligent  materials  and  structures. 

Shape  Memory  Alloys  (SMA)  are  among  candidate  materials  suitable  for  wide  applica¬ 
tion  in  smart  structures.  They  have  the  unusual  material  property  of  being  able  to  sustain 
and  recover  large  strains  (of  the  order  of  10  %  )  without  inducing  irreversible  plastic  defor¬ 
mation  and  to  ’’remember”  a  previous  configuration  and  return  to  it  with  a  temperature 
change.  These  interesting  materieil  characteristics  arise  due  to  distinctive  internal  crys¬ 
talline  transformations  with  temperature  and  applied  stress  (Delaey  et  al.  [1974];  Perkins 
et  al.  [1976];  Funeikubo  [1987];  Wayman  and  Duerig  [1990]). 

The  growing  globjil  interest  in  smart  materied  eind  smart  structure  technology  in  partic¬ 
ular  has  prompted  an  increasing  number  of  investigations  of  SMAs  in  the  past  decade.  The 
result  of  the  research  has  been  increasingly  detailed  information  regarding  the  ciy'stalline 
structure  of  SMA  materials,  a  greater  understanding  of  macroscopic  SMA  material  behav¬ 
ior,  the  development  of  new  alloys  and  processing  techniques,  and  a  dramatic  increase  in 
the  number  of  applications  studied,  which  now  span  a  wide  range  of  products  and  devices. 

The  most  established  commercial  application  of  SMAs,  that  of  connectors  for  hydraulic 
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tubing  in  aircraft,  is  currently  being  overshadowed  by  the  incorporation  of  SMAs  into 
critical  roles  in  a  large  array  of  applications  including  active  vibration  control  of  structures, 
heat  engines,  orthodontic  wire  and  automatic  switches  in  homehold  appliances  (Banks 
Weres  [1976];  Funakubo  [1987];  Rogers  et  al.  [1989];  Falcioni  [1992]). 

Since  Nitinol’s  discovery  in  the  late  1950’s,  a  great  deal  of  research  has  been  reported  on 
the  behavior  of  shape  memory  alloys,  possible  application  and  the  physics  associated  with 
shape  memory  effect.  Rodgers  Robershaw  [1988]  have  demonstrated  that  SMA  materials 
can  be  embedded  in  composite  structures  to  actively  control  and  modify  the  dynamic 
and  structural  behavior  of  SMA  hybrid  composite  lannnates.  Waym2in  [1990,1992,1993] 
discusses  various  practical  applications.  There  have  been  severzJ  attempts  to  construct 
phenomenological  models  capable  of  representing  the  SMA  macrobehavior.  Representative 
works  are  by;  Abeyaratne  and  Knowles  [1993],  Brandon  and  Rogers  [1992],  Brinson  [1993], 
Cory  and  McNichols  Jr.  [1985],  Falk  and  Konopka  [1990],  Ivshin  and  Pence  [1993a,  1993b], 
Liang  and  Rogers  [1990,1992],  Muller  amd  Xu  [1991],  Patoor,  Eberhardt  and  Berveiller 
[1988],  Tanaka  et  al.  [1082,1985,1986,1992],  Tobushi  et  al.  [1991],  Wilmanski  [1993]. 

In  a  recent  review  article  Huo  and  Muller  [1993]  discussed  the  nonequilibrium  thermo¬ 
dynamics  of  pseudoelasticity.  Finite  element  amalysis  of  SMA  behavior  for  one-dimensional 
applications  wais  discussed  by  Brinson  &  hammering  [1993]. 

One  advantage  of  shape  memory  adloys  over  other  types  of  mechanical  or  electrical 
devices  is  that  their  physicail  configuration  can  be  easily,  precisely  and  repeatedly  controlled 
by  often  small  temperature  changes.  In  some  applications,  where  temperature  change  is  in 
itself  the  motivation  for  movement  of  a  mechanical  device,  the  SMAs  can  be  designed  such 
that  no  power  source  is  needed  to  activate  their  motion.  In  other  cases,  a  single  or  series 
of  SMA  wires  can  easily  and  inexpensively  be  given  a  small  temperature  change  which 
initiates  a  phase  transformation  and  consequently  results  in  the  desired  motion  and/or 
stress.  Their  ability  to  achieve  large  strains  near  instantaneously  enables  the  design  of 
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structures  capable  of  extremely  large,  recoverable  deflections.  In  addition,  shape  memory 
alloys  are  relatively  lightweight,  biocompatible,  easy  to  manufacture  and  have  a  high  force 
to  weight  ratio  (Wayman  [1980]). 

Given  the  variety  of  potential  uses  for  SMAs  and  the  high  interest  in  developing  new 
applications,  the  ability  to  accurately  model  and  analyze  structures  containing  SMA  com¬ 
ponents  via  a  finite  element  procedure  is  extremely  attractive,  and  we  Eire  unaware  of  any 
existing  numerical  package  capable  of  addressing  3-D  thermomechanical  behavior  of  shape 
memory  alloys.  The  incorporation  of  the  SMA  numerical  simulation  procedure  into  the 
design  stages  of  new  products  could  reduce  development  times  Eind  costs  dramatically. 
Since  the  properties  of  a  particular  alloy  can  be  easily  and  drasticadly  altered  in  the  man¬ 
ufacturing  process,  the  properties  of  the  SMA  component  in  a  given  design  can  be  varied 
systematically  in  an  automated  analysis  before  production.  This  optimization  procedure 
will  enable  use  of  shape  memory  alloy  components  with  specificEdly  tailored  properties  that 
will  realize  their  full  potential  in  each  individual  application. 

1.2  A  brief  account  of  thermomechanical  behavior  of  Shape  Memory  Alloys 

(Huo,  Muller  [1993]) 

Shape  memory  alloys  in  uniaxial  tensile  tests  exhibit  a  strongly  temperature-dependent 
behavior.  At  high  temperatures  a  typical  load-deformation  curve  is  schematically  shown 
in  Figure  1.1.  The  behavior  characterized  by  this  (P,d)  curve  is  cadled  pseudoelasticity.  It 
is  elastic,  because  the  tensile  specimen  returns  to  the  origin  after  unloading,  but  it  is  only 
"pseudo”  elastic,  because  there  is  a  hysteresis  loop. 

A  hysteresis  loop  indicates  energy  dissipation  or  entropy  pro  hiction,  in  the  loading¬ 
unloading  process,  and  therefore  —  from  the  thermodyneimic  point  of  view  —  the  process 
is  irreversible.  In  other  words,  during  the  process  the  specimen  does  not  remain  in  equi¬ 
librium;  rather  it  runs  through  non-equilibrium  states  for  part  of  the  process. 

In  the  well- developed  non-equilibrium  thermodynamic  theory,  four  "mechanisms”  of 
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Fig.  1.1  (P,d)  diagram  of  shape  memory  alloy  (schematic  picture) 

energy  dissipation  are  commonly  investigated:  heat  conduction,  viscosity,  diffusion  and 
chemical  reaction.  None  of  these  is  present  here  to  a  significant  degree,  but  there  is  a 
phase  transition,  the  so-called  martenistic  transformation.  In  the  unloaded  body  at  high 
temperature  the  crystal  structure  is  cubic  and  the  phase  is  called  austenite  (A).  The  low 
temperature  phaise  has  less  symmetry,  and  it  is  called  martensite  (M).  A  load  can  force  a 
specimen  to  become  martensitic  even  at  high  temperatures  and  Figure  1.1  demonstrates 
this  phenomenon:  The  specimen  is  in  the  A-pheise  on  the  elastic  branch  through  the  origin 
amd  in  the  M-phase  on  the  straight  steep  branch  on  the  right.  The  energj’  dissipation, 
or  entropy  production,  indicated  by  the  hysteresis  occurs  during  the  transition  from  one 
phase  to  the  other. 

Most  of  the  works  on  phase  transitions  have  been  concerned  with  equilibrium  conditions 
of  the  phase  mixture.  Therefore  hysteresis  is  not  considered. 

Close  observation  of  the  phaise  transition  indicates  that  a  specimen  consisting  of  a 
mixture  of  phases  is  composed  of  several  regions  of  pure  phase.  Figure  1.2  shows  a  typical 
schematic  picture  of  a  single  crystal  sample  of  shape  memory  alloy  in  a  tensile  test.  During 
the  transition  from  one  phase  to  the  other,  more  and  more  regions  of  the  new  phase  are 
being  created  and  grow.  At  the  same  time  the  regions  of  the  old  phase  decrease  in  number 
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Fig.  1.2  Single  crystal  sample  of  shape  memory  alloy  composed  of  regions 
of  different  phases  (schematic  picture) 

and  size  and  eventually  they  disappear  completely.  Therefore  the  transition  proceeds  by 
nucleation  and  growth  of  the  new  phase  within  the  old  one.  Since  the  crystal  structures 
and  the  macroscopic  shapes  of  the  two  phases  are  different,  adjacent  regions  of  two  phases 
do  not  fit  each  other  well.  Therefore  additional  energy  is  required  for  the  coexistence  of 
phases,  similar  to  the  surfzu:e  tension  of  a  droplet  surrounded  by  its  vapor. 

In  the  stage  of  nucleation,  this  energy  may  be  calculated  approximately  by  the  use 
of  Eshelby’s  solution  of  inclusions  [1957],  if  the  shape  of  the  nucleus  is  known,  see  e.g. 
the  works  of  Ling  &  Owen  [1981]  and  of  Deng  and  Ansell  [1990]  for  the  case  of  a  shape 
memory  alloy  and  the  work  of  Liu  [1992]  for  a  general  discussion.  Grujicic  [1983],  Olson 
&:  Cohen  [1986]  studied  the  motion  of  a  single  interface  in  a  load  controlled  experiment 
and  determined  the  influence  of  the  defects  at  the  interface  and  inside  the  specimen  on  the 
motion  of  such  an  interface. 

However,  during  the  transition  process  the  places  of  nucleation  are  many,  typically 
of  order  10*,  and  they  are  rather  randomly  distributed  in  the  seimple.  Therefore  it  is 
extremely  difficult  to  obtmn  relevant  results  for  the  macroscopic  behavior  directly  from 
the  study  of  a  single  nucleus  or  interface  emd  approximations  are  unavoidable. 

Load- deformation  diagrams  of  pseudoelasticity  reported  in  the  literature  are  quite  dif¬ 
ferent,  even  for  single  crystals  of  the  seune  composition.  Different  authors  obteiin  very 


diiferent  results,  possibly  because  of  dfferen  t  methods  of  prepeiration  of  the  specimen  es¬ 
pecially  concerning  the  heat  treatment.  Therefore  it  seems  that  at  the  present  state,  it  is 
not  realistic  to  expect  a  good  model  which  can  cover  all  experimental  results. 

For  an  ideal  pseudoelastic  behavior  the  elastic  branches  of  the  two  pure  pheises  are 
parallel  linezu'  functions  of  the  deformation  and  the  hysteresis  loop  is  composed  of  horizontal 
yield  and  recovery  lines.  Figure  1.3  shows  a  typicail  experimental  plot  of  hysteresis  and  its 
idealized  form. 

We  believe  that  ideal  pseudoelasticity  should  be  studied  first,  and  that  its  study  should 
help  us  to  understand  the  essential  features  of  pseudoelasticity.  A  similar  situation  heis 
occurred  in  the  field  of  plasticity;  the  linearly  elastic,  perfectly  plastic  material  has  the 
simplest  plastic  behavior  and  its  study  was  the  first  step  for  a  satisfactory  description  of 
plasticity. 


Strain  [•/.] 

Fig.  1.3  Pseudoelasticity,  left:  experimental,  right:  ideal 

1.2  The  Phenomena  in  Ideal  Pseudoelasticity 

In  order  to  fix  ideas  about  the  problems  we  want  to  investigate,  we  give  a  brief  summary 
of  the  experimental  results  obt2iined  by  Fu  and  Xu  [1992]. 
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A  SMA  single  crystal  specimen  (Cu  Zn  At)  undergoes  phase  transitions  from  austenite 
to  martensite  by  cooling  and  frt)m  martensite  to  austenite  by  heating.  Let  A /  (austenite 
finish)  denote  the  temperature  above  which  the  unloaded  body  is  in  the  A-phase.  An 
isothermal  unizocial  deformation-controlled  tensile  test  at  T  i  Af  exhibits  ideal  pseudoe¬ 
lastic  behavior  in  (P,d),  (d,T)  and  (P,T)  diagrams.  We  proceed  to  discuss  all  three  of 
them. 

First  the  (P,d)  diagram:  In  this  diagram  the  loading- unloading  curve  consists  of  two 
pauallel,  linearly  elastic  branches  and  horizontal  yield  and  recovery'  lines  in-between.  If  we 
interrupt  the  yield  before  it  arrives  at  the  second  elastic  bramch  and  proceed  to  decrease 
deformation,  the  state  of  the  specimen  first  moves  down  steeply  into  the  hysteresis  loop 
along  a  line  parallel  to  the  two  elastic  brzmches,  amd  then  begins  to  recover  at  a  constant 
load  as  shown  in  Figures  1.4. a  and  1.4.c.  Similarly  Figures  1.4. b  and  1.4. d  show  the  internal 
yield  after  the  recovery  process  wars  interrupted.  It  seems  that  there  exists  a  diagonal  line 
inside  the  hysteresis  loop,  as  indicated  by  the  dashed  lines  in  Figures  1.4.c  and  1.4.d  which 
determines  the  turning  points  of  the  internal  loading-unloading  process.  We  shall  call  this 
line  the  diagonal. 

If  we  interrupt  the  yield  or  the  internal  yield  and  unload,  but  reloeid  before  we  arrive  at 
the  diagoneil,  the  reloading  line  is  identical  to  the  unloading  line.  The  specimen  will  resume 
the  yield  however  when  it  returns  to  the  load  where  the  yield  weis  interrupted  before,  as 
if  it  had  remembered  that  load.  A  similar  behavior  is  recorded  when  we  interrupt  the 
recovery  process.  We  cedi  the  steep  lines  inside  the  hysteresis  loop  the  internal  elastic  lines 
and  the  fact  that  the  specimen  yields  or  recovers  only  at  the  load  at  which  it  previously 
crossed  the  diagonal  we  shall  refer  to  as  history  dependence.  Figure  1.5  illustrates  this 
complex  behavior  inside  the  hysteresis  loop. 

When  the  temperature  is  increased  both  the  yield  and  recovery  loads,  Py  and  Pr. 
respectively,  increaise  linearly  with  temperature.  However,  the  difference,  Py  —  Pr,  and 
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Fig.  1.4  Internal  recovery  and  yield.  The  diagonal.  Top:  experimental,  bottom:  ideal 

the  transformation  strain  which  is  the  horizontal  distance  of  the  two  eleistic  branches, 
remain  constzmt.  Thus  the  size  of  the  hysteresis,  i.e.,  (P^  —  Pr)^i,  is  independent  of  the 
temperature,  see  Figure  1.6. 

Figure  1.7  illustrates  how  the  history  dependence  is  influenced  by  a  change  of  temper¬ 
ature.  The  state  of  the  specimen  has  first  arrived  at  the  point  Xz  by  yield  and  elastic 
unloading  at  temperature  Tj,  through  Xj  —  >  X2  —  >  X3.  Then  the  specimen  is  heated 
to  T2  at  fixed  deformation  and  it  remains  at  X3.  Now  upon  increasing  the  deformation 
the  state  moves  first  elastically  to  X4  and  then  begins  to  yield  at  a  load  higher  than  the 
previous  one.  The  dashed  lines  indicate  the  yield  and  recovery  lines  of  the  body  at  T2. 

Second  the  (d,T)  diagram:  This  diagram  is  taken  by  chainging  the  temperature 
of  the  specimen  which  is  under  a  fixed  load.  We  neglect  thermal  expansion,  so  that  the 
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Fig.  1.5  Internal  elzisticity  and  history  dependence.  Top:  experimental,  bottom:  ideal. 

(d,T)  diagram  exhibits  two  horizontal  Unes,  of  which  the  higher  one  is  for  the  M- phase 
and  the  lower  one  for  the  A-phase,  and  two  vertical  jiimps  where  the  phase  transitions 
occur.  Similar  to  the  (P,d)  diagram,  internal  jumps  are  observed  as  shown  in  Figures  1.8. a 
and  1.8.C,  and  they  determine  a  line,  the  dashed  line  in  the  figure,  inside  the  hysteresis 
loop,  which  corresponds  to  the  diagonal  in  Figure  1.4.  Figures  1.8.b  and  1.8.d  show  the 
thermal  path  which  is  obtained  by  interrupting  the  jump  from  the  A-phase  to  the  M-phase, 
then  heating  and  re-cooling  the  sample  before  it  arrives  at  the  diagonal.  The  state  of  the 
specimen  moves  back  and  forth  along  the  same  horizontal  path. 
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Fig.  1.6 


Fig.  1.7 
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Temperature  dependence  of  the  hysteresis,  left:  experimental,  right:  ideal 


strain  1%) 


History  dependence  of  a  temperature  change,  left:  experiment,  right:  ideal 
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Fig.  1.8  Internal  jumps  and  thermal  path.  Top:  experiment,  bottom;  ideal 

Third,  the  (P,t)  diagram:  Here  we  keep  d  constant  and  monitor  the  load  as  the 
temperature  is  changed.  Again,  ignoring  the  tendency  for  thermal  expansion  we  obtain  two 
horizontal  lines;  the  lower  one  is  for  the  M-phase,  while  the  higher  one  is  for  the  A-phase. 
The  transition  from  M  to  A  with  increasing  temperature  is  a  straight  line  inclined  to  the 
T-axis,  and  there  is  a  higher  parallel  line  for  the  reverse  transition,  see  Figure  1.9. 

By  interrupting  the  A  —  >  M  transition  and  increasing  temperatures  the  state  of  the 
specimen  moves  horizontally  into  the  hysteresis  loop  on  a  reversible  path  as  indicated  in 
Figures  1.9.b  and  1.9.d.  When  the  temperatiire  is  lowered,  interrupting  the  M  —  >  A 
tremsition,  and  the  lowering  proceeds  fax  enough,  there  is  again  evidence  of  a  line  dashed 


in  Figures  1.9.c  said  1.9.d,  which  corresponds  to  the  diagonal  of  Figure  1.4.  On  that  line 
the  specimen  begins  the  phase  transition  inside  the  hysteresis  loop  and  pwallel  to  the  outer 
trsmsition  lines  sis  shown  in  Figures  1.9. a  and  1.9.c. 

The  word  pseudoelasticity  often  refers  only  to  the  behavior  in  the  (P,d)-diagrams  ob- 
tsuned  by  T  —  >  Aj.  However,  if  the  load  is  reasonably  big,  the  M-phase  obtained  in 
a  (d,T)  experiment  is  of  only  one  variant,  said  we  may  treat  (d,T)-diagrams  and  (P,T)- 
diagrams  by  using  the  same  approach  as  for  the  (P,d)-diagrams.  We  may  even  construct 
the  (d,T)-  and  (P,t)-diagrams  directly  from  a  series  of  (P,d)-diagrsims  obtained  at  different 
temperatures. 
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Fig.  1.9  The  (P,T)-diagram.  Top:  experiment,  bottom:  ideal. 
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CHAPTER  2 


A  CONSTITUTIVE  MODEL  OF  SHAPE  MEMORY  ALLOYS: 
PSEUDOELASTIC  BEHAVIOR 
AND  ITS  ALGORITHMIC  IMPLEMENTATION 

2.1  Introduction 

The  purpose  of  the  present  development  is  to  model  the  so-called  pseudo-eljistic  behav¬ 
ior,  represented  in  Figures  1.3  and  1.4  which  occurs  in  shape  memory  alloys  (SMA)  under 
stress  loading-unloading  pattern  at  constant  temperature.  In  particular  we  will  present  a 
constitutive  model  able  to  predict  the  stress-induced  phase  treinsformation  &om  austenite 
to  martensite  during  loading  and  the  reverse  transformation  from  martensite  to  austenite 
under  unloading.  We  are  also  concerned  with  an  appropriate  description  of  the  mate¬ 
rial  behavior  when  incomplete  transformations  occur.  These  incomplete  transformations 
usually  induce  internal  loops,  as  qualitatively  represented  in  Figures  1.4  and  1.5. 

Since  the  present  discussion  is  limited  only  to  the  pseudo-eleistic  behavior  of  shape 
memory  Jilloys,  we  consider  only  iso- thermal  processes. 

The  present  section  is  organized  as  follow:  we  start  discussing  the  basic  assrunptions  and 
describing  the  continuous  model,  with  the  appropriate  phase  transitions.  We  present  the 
corresponding  discrete  time  model  and  discuss  the  relative  algorithmic  implementation, 
based  on  the  use  of  a  radial  return  map  scheme  for  the  integration  of  the  constitutive 
equations.  We  also  present  the  form  the  algorithmic  tangent  consistent  with  the  discrete 
model.  We  conclude  the  paper  with  some  numerical  examples. 

2.2  Basic  Assumptions  and  Continuous  Time  Model 

Limiting  the  discussion  to  the  realm  of  small  deformations,  we  may  split  the  toted  strain 
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(2.1) 


e  into  an  elastic  component  c'*  ,  and  an  inelastic  component  e': 

e  =  c*'  + 

The  inelastic  component  c’  is  assumed  to  be  induced  by  possible  phase  transformations 
the  material  may  imdergo.  To  simplify  the  notation,  the  dependency  of  the  variables  on 

the  time  t  is  not  explicitly  stated.  As  usual,  we  may  also  split  the  stress  O'  and  the  strain 

e  in  the  volumetric  eind  the  deviatoric  components; 

<r  =  »  pi  (2.2) 

e  =  e  +  ^ei  (2.3) 

where:  p  =  jtr(<r)l  is  the  pressure  and  9  =  tr(c)  is  the  volume  change.  In  particu¬ 
lar,  recalling  equation  2.1  we  may  split  both  the  elastic  and  the  inelastic  strmn  into  the 
volumetric  and  the  deviatoric  components; 

e  =  cO  +  +  1  (2.4) 

Basic  assumptions  under  which  we  develop  the  model  are: 

•  the  phase  transformations  depend  only  on  the  deviatoric  part  of  the  stress  and  affect 
only  the  deviatoric  part  of  the  strain, 

•  in  the  absence  of  phase  transformations  the  material  has  a  linear  elastic  behavior. 
Accordingly,  9’  =  0  auid  the  elastic  relations  may  be  written  eis: 

p  =  KB  (2.5) 

s  =  2Ge''  =  2G(e  -  e‘)  (2.6) 

where  K  is  the  bulk  modulus  and  G  is  the  shear  modulus.  To  be  consistent  with  the 

notation  usually  adopted  in  the  development  of  uniaxial  model  for  SMA  (Liang  k.  Rogers 
[1990],  Brinson[1993]),  we  rescale  the  inelastic  contribution: 

c’  =  (2.7) 
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where  is  a  material  parameter,  measurable  experimentally  and  related  to  the  so-called 
maximum  residual  strain  under  uniaxial  state  ei  through  the  relation: 

e'i  =  yfet  (2.8) 

The  evolutioneiry  equation  for  *  is  assumed  in  the  form: 

X  =  in  (2.9) 


where: 


n  = 


(2.10) 


is  a  unit  vector  in  the  stress  direction.  In  plasticity  literature,  n  is  a  unit  vector  perpen¬ 
dicular  to  the  yield  function.  ^  is  a  scalar  internal  variable  representing  the  martensite 
fraction.  A  super-posed  dot  indicates  time  derivative.  Due  to  the  interpretation  of  we 
require  that  0  ^  <  1  and 


{ 


^  =  0  •O  the  material  is  all  austenite, 

^  =  1  the  material  is  all  mEU'tensite, 


(2.11) 


We  are  now  left  with  the  need  of  describing  the  kinetics  of  the  phase  transformations,  i.e. 
we  have  to  specify  the  rules  for  the  activation  euid  evolution  of  the  pheise  transformations. 
Clearly,  we  have  two  different  transformations  to  consider: 

1.  production  of  martensite,  which  means  conversion  of  austenite  into 
martensite  {A  — »  M), 

2.  production  of  austenite,  which  means  conversion  of  martensite  into 
austenite  (AI  — ►  A). 

It  h^ls  been  shown  experimentally  that  under  uniaxial  loading  conditions,  each  phase  trans¬ 
formation  occurs  within  an  appropriate  range  of  stresses;  in  particular,  indicating  with  a 
the  uniaxial  stress,  we  have: 


,  <7’>0  A  —*  M 

<^Aj  <  <  <^A,^  d  <  0  =>  M  —*  A 
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where  it  is  clearly  stated  that  the  conversion  of  austenite  into  martensite  (A  —*  M)  occurs 
during  loading,  while  the  conversion  of  martensite  into  austenite  (M  —*  A)  occurs  under 
unloading. 

Recalling  that  the  transformations  are  assumed  to  be  only  function  of  the  deviatoric 
part  of  the  stress,  we  may  generalize  the  concept  of  unieixial  stress  ramge  using  a  J2-type 
norm  where  J2  is  the  second  invariant  of  the  deviatoric  part  of  the  stress.  In  particular, 
we  set: 

S  =  M  =  (2.13) 


and  as  conditions  for  the  phase  transformations  we  require: 

Sm,  ^  S  >  Q  a  —*  M 

Sa,<S<Sa.,  S  <0  M  A 


(2.14) 


The  three- dimensioned  S-bounds  are  related  to  the  uniaxial  (7-bounds  through  the  relations: 

Sm,  =  ,  Sm.  = 


(2.15) 


2.3  Austenite-Martensite  Phase  Transformation 

We  eissume  a  linear  variation  of  the  martensite  fraction  (  from  0  to  1  within  the  stress 
rEinge  [5a/,,  vmder  loeuiing  conditions.  Accordingly,  the  following  equation  must  be 

identically  satisfied: 

1  -  C’’ 

^am[0  =  i  - 

where: 

^AM  —  Smi  -I-  5Af, 

Ham  =  Smj  —  Sm. 

and  is  the  maximum  value  obtained  by  the  martensite  fraction  during  the  previous 
(reverse)  martensite-austenite  transformation. 
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Ham 


1  +  P 


(2.16) 
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2.4  Martensite-Austenite  Phase  TVansformation 


We  assume  a  linear  variation  of  the  martensite  fraction  ^  from  1  to  0  within  the  stress 
range  [S^,,  S^J,  under  xmloading  conditions.  Accordingly,  the  following  equation  must  be 
identically  satisfied: 


^Ma{0  = 


2 


—25  + 
Hma 


+  1 


=  0 


where: 


f^MA  =  Sa,  +  Sa, 


(2.17) 


M  A  —  ^Af  ^A, 

and  is  the  maximum  value  obtained  by  the  martensite  fraction  during  the  previous 
austenite-martensite  transformation: 


2.5  Internal  Loop  Descriptions 


In  order  to  have  a  correct  description  of  the  internal  loops  for  the  case  in  which  incom¬ 
plete  transformations  occur,  we  modify  the  range  of  the  transformation  depending  on  the 
value  of  the  martensite  fraction  during  the  previous  transformation. 

For  example  (referring  to  Figxires  2.5.1  and  2.5.2)  for  the  case  of  conversion  of  austenite 
into  martensite,  we  set  the  range  within  which  the  tremsformation  occurs  as: 


Vm.  =  CSa.  + 


(2.18) 


Vm,  =  ^Sa.  +  (1-^'’)5m, 

where  is  the  value  attained  by  the  meirtensite  fraction  during  the  previous  reverse 
transformation  (i.e.,  during  the  conversion  of  martensite  into  austenite).  For  the  case  of 
conversion  of  martensite  into  austenite  we  similarly  set: 


va.  =  esA.  +  {i-insM. 


(2.19) 


Va,  =  CSa,  +  (l-r)5Ar. 

where  is  the  value  attained  by  the  martensite  fraction  during  the  previous  reverse 
transformation  (i.e.  during  the  conversion  of  austenite  into  martensite). 


2.6  The  Kuhn- Tucker  Conditions 


Let  functions  ^maH)  be  the  criteria  for  {A  — »  M)  and  {M  — »  A)  phase 

transformations,  respectively,  similar  to  yield  function  in  plasticity.  In  the  rate  independent 
constitutive  relations  for  SMA,  the  requirement  that  ^  <  oo  necessitates: 


:r{0  :=Q  ^ 


(2.20) 


where 


= 


^maU) 
'  :FAM{i) 


^  for  M  — »  A,  (S  >  0) 
<=►  for  A  —*  M,  (5  <  0) 


(2.21) 


and  S  =  ||s||  =  s/TTl  =  s/TJ^-  This  is  the  transformation  criteria  of  rate-independent 
SMA.  The  transformation  criteria  and  the  loading/unloading  conditions  can  be  expressed 
in  Kuhn-Tucker  form  as  follows. 


CASE  1.  Austenite  — ►  Martensite  transformation 


J^amU)  <  0,  ^  >  0;  iff  (7  >  0  and  <7  <  cr*  (2.22) 

where 

a*  —  (TMj  (complete  transformation) 

<7*  =  Vm;  (incomplete  transformation) 

<7jvf;  and  Vmi  are  critical  stresses  for  austenite  to  martensite  trEinsformations  corresponding 
to  outer  and  inner  hysteresis  loops  respectively. 

CASE  2.  Martensite  — *  Austenite  transformation 


^maH)  >  0,  ^  <  0;  iff  <7  <  0  and  cr  >  a  (2.23) 

where 

a  =  a Af  (complete  transformation) 
a  =  V A,  (incomplete  transformation) 


and  Vaj  are  critical  stresses  for  martensite  to  austenite  transformations  corresponding 
to  outer  and  inner  hysteresis  loops  respectively. 

Combining  the  above  expressions,  we  have 


|.F*(oi >0  ^  i  =  o 


(2.24) 


and 

1^1  >  0  =>  =  0 


(2.25) 


These  two  conditions  (as  in  plasticity)  express  the  physiced  requirements  that  the  stress 
must  be  admissible  and  the  phase  transformation,  in  the  sense  of  a  non-zero  martensite 
fraction  rate  ^  0  can  take  place  only  on  the  phase  transformation  surfaces 

Within  the  context  of  plasticity,  these  conditions  are  classical  in  the  convex  mathemat¬ 
ical  programming  literature  and  go  by  the  name  of  Kuhn- Tucker  conditions.  These  two 
equations,  (2.24,  2.25)  result  in 

r{i)i  =  o  (2.26) 

In  computational  plaisticity  literature,  this  equation  is  referred  to  as  the  “consistency  con¬ 
dition.” 

We  can  define  the  phase  transformation  surface  (known  as  the  yield  surface  in  plasticity) 
as 

dJ^\i{s))  =  {<T€  ^|:F*(e(s))  =  0}  (2.27) 

For  the  case  when  the  state  lies  on  the  phase  transformation  surface 


and  set 


1^1  >0  if 


only  if  .r*(^(s))  =  0 

(2.28) 

^am{0 

<  0 

or 

(2.29) 

^maH) 

>  0 

23 

BASE 


We  have  an  additional  condition. 

P{()i  =  0  (2.30) 

This  condition  is  referred  to  as  the  “persistency  condition,”  and  corresponds  to  the  physical 
requirement  that  in  order  for  (  to  be  non-zero  (i.e.,  |^|  >  0),  the  stress  state  a  € 
must  ’’persist”  on  d^*(((s))  so  that  ^*{0  =  0. 

2.7  Discrete  Time  Model  and  Algorithmic  Implementation 

We  now  present  the  discrete  time  counterpau’t  for  the  model  discussed  above,  paying 
particular  attention  to  an  implementation  within  a  radial  return  map  algorithm.  The  form 
of  the  tangent  tensor  consistent  with  the  discrete  model  is  also  addressed. 

2.7.1  Discrete  Equations  and  Integration  Algorithm 

FVom  a  computational  standpoint  we  treat  the  non-linear  behavior  of  a  material  as 
a  strain  driven  problem,  since  in  a  finite  element  implementation  the  stress  history  is 
computed  from  the  strain  history  by  an  integration  technique,  such  as  a  return  mapping 
algorithm.  Accordingly,  we  now  introduce  a  discrete  counterpart  of  the  equations  presented 
ezu-lier  and  review  the  integration  algorithm. 

Let  [0,  T]  C  ^  be  the  time  interval  of  interest  and  consider  two  time  values  within  it, 
say  tn  and  <„+i  >  t„,  such  that  is  the  first  time  vailue  of  interest  after  t„.  To  minimize 
the  appearance  of  subscripts  (aj'd  to  make  the  equations  more  readable),  we  introduce  the 
convention: 

On  =  a  =  a(t„+,) 

where  a  is  any  generic  quantity.  Accordingly,  in  the  discrete  time  setting  the  subscript  n 
indicates  a  quantity  evaluated  at  time  t„,  while  no  subscript  indicates  a  quantity  evaluated 
at  time  tn-i-i  • 

We  assume  that  the  solution  is  known  at  time  t„  and  given  by  the  state; 


We  wish  to  compute  the  solution  at  time  t„+i,  given  the  strsdn  c.  Using  a  backward  Euler 
integration  formula  for  the  scaled  inelastic  strain,  we  obtain; 

*n+i=*n  +  An  (2.32) 

where 


A  = 


=►  ^  +  A 


(2.33) 


Substitution  of  equations  (2.29)  and  (2.30)  into  equation  (2.5)  yields: 

8  =  2G[e  -  e\^x„]  -  2GelXn  (2.34) 


In  the  above,  A  is  an  unknown  quantity  and  is  computed  by  means  of  an  integration 
algorithm,  such  as  a  retvu-n  mapping.  Initially  suggested  by  Maenchen  and  Sack  [1964], 
the  return  mapping  algorithm  provides  an  efficient  and  robust  integration  scheme,  based 
on  a  discrete  enforcement  of  the  limit  equation.  It  belongs  to  the  fEimily  of  elastic-predictor 
inelastic-corrector  algorithms  and,  hence,  is  a  two  part  algorithm.  In  the  first  part,  a  pmely 
elastic  trial  state  is  computed;  in  the  second,  if  the  trial  state  violates  the  materiail  model 
constitutive  equation,  a  correction  is  computed  using  the  trial  state  as  initiaJ  condition  and 
applied  such  that  the  finzJ  state  is  fully  consistent  with  the  discrete  model.  The  algorithm 
has  been  widely  studied  (Nagtegal  [1982],  Simo  &  Taylor  [1985],  Simo  &  Hughes  [1993])  as 
has  its  stability  (Krieg  &  Krieg  [1977],  Simo  &  Govindjee  [1991]).  Additional  discussion 
of  the  algorithm  and  its  theoretical  implication  can  be  found  in  the  literature. 

We  shall  now  discuss  the  two  steps  of  the  algorithm  in  more  details. 


•  IVial  state:  we  assume  that  in  the  interval  [<„,  <n+i]  no  phase  transformation  occurs 
(i.e.  X  =  Xni  which  implies:  A  =  0).  As  a  result,  we  have  as  trial  values: 


=  0 

(2.35) 

=  »n 

(2.36) 

=  2G  [e  - 

(2.37) 
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If  the  elastic  trial  state  is  admissible,  i.e.  it  does  not  satisfy  the  conditions  for  a  phase 
transformation,  then  it  represents  the  new  solution  at  f„4.i  and  the  second  part  of 
the  algorithm  is  skipped.  If  the  elastic  trial  state  is  not  admissible,  a  correction 
hsis  to  be  performed,  i.e.  the  phase  transformation  must  be  taken  into  account  to 
determine  the  real  solution  state. 

•  Inelastic  correction:  enforcing  the  satisfaction  of  the  phzise  trzmsformation  model, 
the  parameter  A  may  be  computed,  as  shown  for  both  transformations  {A  M 
and  M  —*  A)  later.  Equations  (2.32)  and  (2.34)  can  be  now  rewritten  in  terms  of 
the  trial  state  and  A: 

*  =  a;*’’  +  A  n  (2.38) 

s  =  s‘’'-2Ge‘jr  An  (2.39) 

which  allow  us  to  compute  the  inelastic  solution.  Since  the  phase  transformations 
occur  within  specific  upper  and  lower  bound,  before  updating  the  solution  with  the 
inelastic  one  described  in  equations  (2.31)  and  (2.32),  it  is  important  to  test  that 
such  solution  is  still  in  the  range  of  the  pheise  transformation. 

For  both  transformations  considered  in  this  work,  the  parameter  A  may  be  computed 
solving  only  a  scalar  equation.  In  fact,  we  observe  that  s  is  by  definition  in  the  n-direction, 
i.e.,  8  =  l|«||n.  Consequently,  a  scalar  relation  between  the  norms  of  8  and  «*'’  may  be 
generated; 

S=  ||s||  =  ||*‘n|-2Ge'iA  =  5‘"-2Ge'iA  (2.40) 

2.7.2  Discrete  Phase  Transformations 

Substitution  of  equations  2.30  and  2.34  into  the  constitutive  equations  for  both  the 
transformations  returns  the  corresponding  discrete  models  in  the  form; 

=  U  +  -  Ram  -  2Ge‘^A)  +  =  0  (2.41) 

^MA  =  +  A  —  RMAi^S^^  —  2Ge^A)  +  RmaSa,  =0  (2.42) 
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where: 


1  _  £P  fP 

Ram  =  -7r - 1  Rma  =  -  tt —  (2-43) 

^AM  tlMA 

Due  to  the  linearity  of  such  equations,  the  parameter  A  may  be  computed  in  closed  form 
for  both  the  transformations  and  it  is  given  by: 


A-^M  ^ 
M  —*  A 


RamS*’^  +  -  £„ 

1  +  2  G  Ram 
Kma{S^^  -  Sa,)  -  in 
1  +  2Gc^  Rma 


(2.44) 

(2.45) 


2.7.3  Conditions  for  the  Activations  of  the  Phase  TVansformations  in  the  Dis¬ 
crete  Model 

The  conditions  for  the  activations  of  the  phase  transformations  described  in  Section  2.5 
for  the  continuous  time  model  must  now  be  rewritten  for  the  corresponding  discrete  time 
model.  The  first  test  we  perform  is  a  check  on  the  loading-unloading  condition,  which  in 
a  discrete  model  can  be  rewritten  as: 


(2.46) 


>  Sn  =>  loading  =>  possible  A  M  transformation 
5**^  <  5n  =►  unloading  possible  M  —*  A  transformation 
After  checking  which  transformation  is  possibly  active,  we  compute  the  corresponding 
A  parameter  (accordingly  to  equations  2.38  or  2.39)  and  check  that  the  final  inelastic 
candidate  state  still  satisfied,  it  means  that  the  final  state  is  still  in  a  phase  treinsition 
range  and  we  can  update  the  solution. 


2.7.4  Discrete  Tangent  Tensor 

We  address  the  form  of  the  tangent  tensor,  consistent  with  the  discrete  model.  The  use 
of  a  consistent  tangent  tensor  preserves  the  quadratic  convergence  of  a  Newton  method, 
which  we  adopt  for  the  incremental  solution  of  a  finite  element  scheme.  We  start  from  the 


1 

1 

linear  elastic  relation  between  8  and  e: 

1 

8  =  2G  [e  -  cj^asn]  -  2Gt\Xn 

(2.47) 

1 

and  by  linearization  we  obtain; 

1 

d8  =  2Gde  —  2Ge*idXn  —  2G  e*i\dn 

(2.48) 

■ 

1 

Keeping  in  mind  that; 

8  8 

"  “  11*11  ■  [(.  :  .)S' 

(2.49) 

1 

we  can  compute  its  variation; 

1 

(2.50) 

■ 

where  the  fourth  order  tensor  N  is  the  orthogonal  projection  operator  on 

the  plaine  with 

■ 

unit  normal  n,  such  that; 

1 

Nn  =  0,  NN  =  N 

(2.51) 

1 

Accordingly,  equations  2.41  can  be  rewritten  as; 

■ 

■ 

[I  +  aN]d8  =  2Gde  -  2Ge[dXn 

(2.52) 

1 

1 

1 

where; 

2Ge\X 

“  "  Ikll 

The  equation  can  be  solved  for  ds: 

■ 

■ 

d8  =  2G[I-\-bN]de  -  2GeldXn 

(2.53) 

1 

1 

where: 

b  =  ° 

1  -1-  a 

1 

Assuming  that  from  the  linearization  of  the  phase  transformation  equation  we  obtain  a 

■ 

relation  of  the  type  dX  =  A\n  :  de]  (as  we  will  prove  in  the  next  section), 

we  C2in  get  the 

1 

1 
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incremental  relation  between  the  total  stress  tr  and  the  total  strain  e,  consistent  with  the 
discrete  model; 

da-  =  Dditc  de  (2.54) 

where  the  algorithmic  tangent  tensor  is  given  by: 

Ddisc  =  K{l(S)l)  +  2G{l-C)Idtv  +  2G{C -A){n®n) 

where 

Idev  =  /  (2.55) 

We  reczJI  that  A  comes  from  the  linearized  limit  equation  for  the  specific  material  model, 
while: 

2££^ 


2.7.5  Linearization  of  the  Phase  Transition  Equations 

We  start  form  the  discrete  equation  for  the  austenite-martensite  phase  transformation: 


=  {„  +  a  -  RamS  +  - 

T.  =  0 

2 

(2.56) 

and  by  linearization  we  get: 

^^AM  =  dX  —  Ram  dS  =  0 

(2.57) 

Since: 

dS  =  2Gn  :  de  —  2GeidX 

(2.58) 

we  get; 


Similarly,  starting  from; 


d\ 


2GRam  j 

l  +  2GelRAM'''  ^ 


^MA  =  +  A  -  RmA  -5  =  0 


(2.59) 


(2.60) 


by  linearization  we  get: 


d\ 


2GRma  , 

\  +  2GelRMA'^'  ® 
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(2.61) 


2.8  Numerical  Examples 

To  test  the  integration  algorithri  for  our  model,  we  run  some  tiniaxial  test  under  dis¬ 
placement  control.  The  material  properties  are: 

E  =  lOkPa,  1/  =  0.3  (2.62) 

while  the  transformation  range  are: 

OM,  =  10,  OM,  =  15,  CA,  =  20,  a  A,  =  25  (2.63) 

with: 

Cm  =  1,  =  1,  Cl  =  3  (2.64) 

The  test  are  organized  as  follow: 

•  Test  1;  we  run  a  loading- imloading  test  with  incomplete  A  —*  M  transformation  and 
complete  M  —*  A  transformation  (see  Figure  2.8.1); 

•  Test  2:  we  run  a  loading-unloading  test  with  complete  A  M  transformation  and 
incomplete  M  —*  A  transformation  (see  Figure  2.8.2); 

•  Test  3;  we  run  a  loading- unloading  test  with  incomplete  A  —*  M  transformation  and 
incomplete  M  —*  A  transformation  (see  Figure  2.8.3); 


T 


e 


e 


Figure  2.S.1  Picudo-elastic  behavior:  modification  of  the 
ff-range  for  the  conversion  of  austenite  into  martenaite. 


Figure  2.5.2  Pteudo-elattic  behavior:  modification  of  the 
«-range  for  the  oonvenion  of  martenaite  into  austenite. 


Fig.  2.8.1  Stress  Vs  strain  plot.  Complete  M  —*  A  and  incomplete  A  M  transformation. 
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CHAPTER  3 


FINITE  DEFORMATION  FORMULATION 
FOR  SHAPE  MEMORY  ALLOYS 

3.1  Introduction 

This  chapter  presents  a  finite  deformation  formulation  for  shape  memory  alloys.  Let 

be  the  deformation  gradient,  where  X  is  the  reference  configuration  and  is  the  mapping 
from  reference  to  current  configuration,  i.e.,  *  =  f(X,t)  =  ft(X).  F  is  assumed  to  be 
s\ifficiently  smooth,  orientation  preserving,  invertible  and  subject  to  the  constraint 

J  =  detF  >  0  (3.2) 

to  ensure  that  the  material  volume  element  remains  positive,  and 

dv  =  detF  dV  =  J  dV  (3.3) 

where  dV  is  a  volume  element  in  the  reference  or  the  material  configuration  and  dv  is  the 
corresponding  form  in  the  current  or  the  spatial  configuration. 

In  the  present  formulation  a  multiplicative  decomposition  of  the  deformation  gradient 
into  eleistic  and  inelastic  parts  is  assumed 


F  =  F'F‘  (3.4) 

where  the  superscripts  e  and  i  represent  elastic  and  inelastic  parts  of  the  deformation 
gradient,  respectively.  The  inelastic  part  is  associated  with  pnase  tremsformation.  See 
Simo  [1988]  for  details.  The  inelastic  part  of  the  deformation  gradient  defines  an  additional 
configuration  Bt,  or  intermediate  configuration.  The  discretized  form  of  (3.4)  is 

e„+,  =  F„V,F'+,  (3.5) 


BASE 


3.2  Phase  Transformation  Flow  Rule 


We  can  derive  the  flow  rule  for  inelastic  or  phase  transformation  strain  by  tedcing  the 
time  derivative  of  the  inelastic  deformation  gradient. 


F'iXj) 


d  (d^'ixj)\ 
dt\  dx  ) 

GRADV’(A',t) 


=  gradt>’(*,  t)  F' 

=  6i,JV(S)f<  (3.6) 


where  V*  is  the  inelastic  material  velocity,  t>’  and  u'  are  the  inelastic  spatial  velocity  and 
displacement  fields,  respectively.  In  the  above  equation  GRAD  V  is  the  ine’^istic  material 
velocity  gradient,  while  gradv*  is  the  inelastic  spatial  velocity  gradient.  Furthermore  we 
have  used  the  finite  strain  coimterpart  of  the  small  deformation  inelastic  strain  (2.7)  to 
arrive  at  eq.  (3.6).  We  can  write  (3.6)  as 

F*F‘“  =  ^eliV(S) 

^InF'  ={ei.N(S) 

Integrating  both  sides  from  to  we  get 


F<n+1  . 

/  dlnF*  =  /  itlN{S)dt 

Jt„ 

InF*^,  -  lnF,j  =  A*ellV(S) 


where  we  have  used  eq.  (2.33)  on  the  right  hand  side  for  A*.  Equation  (3.7)  can  be  written 
as 

InF'+.Fr'  =  Vei,7V(S) 

n+i  F^'  =  exp  (A-  ei  N{S)) 

Consequently, 

F-^,  =  exp(A*e'^N(S))F,;  (3.S) 
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where  A*  is  the  increment  in  the  martensite  fraction  (  and  is  obtained  through  a  return 
mapping  algorithm  via  satisfaction  of  the  corresponding  phase  transformation  criteria,  i.e., 


^ 


^AM  ^  enforcing  =  0 

A"^  =►  enforcing  ^^"*(0  =  0 


The  manner  in  which  the  flow-rule  is  discretized  is  an  essential  part  of  the  method  of 
extension,  and  consists  of  taking  the  inelastic  flow  direction  N  to  be  constant  throughout 
the  increment  and  equal  to  its  final  value  This  reduces  the  phase  transformation 

flow  rule  to  a  system  of  linesu"  equations  for  the  inelastic  deformation  gradient  F'  with 
initial  condition  with  an  exact  solution  given  by  the  exponentiEil  mapping  (3.8). 

Remark:  By  expressing  the  “phase  transformation”  flow  rule  and  the  elastic  response  on 
the  intermediate  configuration,  the  formulation  automatically  satisfies  material  frame 
indifference.  See  Cuitino  et  al.  [1992]. 

3.3  Constitutive  Relations 

Direct  use  of  deformation  gradient  F  complicates  the  development  of  constitutive  equa¬ 
tions  zmd  it  is  common  to  introduce  deformation  measures  which  are  related  completely 
to  either  the  reference  or  the  current  cinfiguration.  Accordingly,  for  the  reference  configu¬ 
ration  the  right  Cauchy-Green  deformation  tensor,  C  is  introduced  as 

C  =  F^F  =  (3.10) 

In  the  reference  configuration,  we  define  the  Green  strain  tensor  E  as 

E  =  \{C-l)  (3.11) 

where  1  is  the  rank  two  identity  tensor  with  respect  to  the  reference  configuration. 

We  can  write  the  elastic  right  Cauchy-Green  transformation  tensor  in  the  intermediate 


configuration  as 


(3.12) 


where  C  is  the  full,  right  Cauchy-Green  deformation  tensor  defined  over  the  intermedi¬ 
ate  configuration.  Accordingly,  we  define  the  Green  strain  tensor  on  the  intermediate 
configuration 

£  =  i(C-l)  (3.13) 

Since 

ilogC  =  i(C  -  1)  -  i(C  -  1)(C  -  1)  +  •  •  •  (3.14) 

So  we  define  another  strain  measure  on  the  intermediate  configuration,  called  the  Henky 
strain 

•®(Henky)  =  2^°®^  (3.15) 

The  second  Piola-Kirchhoff  stress  tensor  corresponding  to  the  Henky  strmn  measure  defined 
over  the  intermediate  configuration  can  be  expressed  as 

■s  =  s  (h°s^)  (316) 

Equations  (3.5),  (3.8)  and  (3.16)  define  a  system  of  nonlinear  equations  which,  for  given 
Fn+i  can  be  solved  for  the  updated  state  variables  and  A.  (Note  A  =  A^.  As 

in  the  small  strain  case,  it  is  possible  to  reduce  the  system  to  a  single  equation  for  A,  which 
can  be  solved  by  using  the  return  mapping  algorithm  embedded  in  the  Newton- Haphson 
iteration  scheme. 

Remark:  The  generalized  strain  measures  intorduced  by  Seth  (1964a,  1964b)  have 
successfully  been  applied  to  the  development  of  material  constitutive  equations  for 
plasticity,  nonlinear  viscoelasticity  and  rubber-like  eleisticity.  In  the  material  descrip¬ 
tion,  these  strain  measures  are  related  to  the  right  Cauchy-Green  deformation  tensor 
C  through  the  expression 

l(C”*/2  _  /)  ^  ^  0 

1  (3.17) 

-logC  m  =  0 
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where  I  is  the  identity  tensor  and  the  exponent  m  may  take  on  noninteger  as  well  as 
integer  values.  For  integer  values  of  m  these  equations  reduce  to  strain  measures  of 
Cauchy  and  Henky  In  fact  the  validity  of  these  equations  (3.17)  for  nonintegreJ  real 
values  of  m  is  assumed  if  —  I)  is  defined  to  be  coaxial  with  the  right  stretch 

tensor  U  =  Euid  having  principal  value  —  1),  (o  =  1,2, 3),  Aq  being  the 

eigenvalues  corresponding  to  U.  As  pointed  out  by  Seth  [1964a,  1964b],  by  introducing 
a  degree  of  freedom  in  the  exponent  m,  it  is  possible  to  condense  the  nonlineau:  effects  of 
deformation  into  the  definition  of  strain,  and  thus  rely  less  on  representing  the  nonlinear 
behaviour  in  the  constitutive  equations. 


3.4  Radial  Return  Algorithm 

Within  the  context  of  radial  return  algorithm,  the  right  Cauchy- Green  deformation 
tensor  for  elastic  predictor,  i.e.,  the  step  in  which  the  phzise  transformation  is  frozen  and 
a  purely  elastic  deformation  is  assumed  is  defined  as 

<?;+!  =  e'r+iF'„+,  {3.18) 


while  is  defined  in  equation  (3.4).  Using  (3.5)  and  (3.8)  with  A  =  0  and  substuting  in 
(3.18) 

Cn+i  F'-’  (3.19) 

This  is  the  large  strain  counterpart  of  the  elastic  predictor  as  can  be  derived  from  (2.37). 
We  can  write  the  Henky  strain  (3.15)  for  the  elastic  trial  state  as 

E{Htnky)„^.t  =  (3.20) 


The  corresponding  second  Piola-Kirchhoff  stress  tensor  defined  on  the  intermediate  con¬ 
figuration  is 


5  = 


(3.21) 
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We  use  a  J2  type  norm  for  the  elastic  trial  state  defined  as 


5*^  =  IlSll  =  v'Sn+i  :  S„+,  (3.22) 

At  this  stage  it  is  checked  if  the  state  is  in  the  phase  transformation  range,  i.e.,  if 

for  all  5>0  (3.23) 

and 

{S^.,Sa,}  for  all  5<0  (3.24) 

then  we  have  the  solution  at  time  t„4i  and  the  second  part  of  the  algorithm  is  skipped. 
However  if 

Sm.  <  <  Sm,  for  5  >  0  (3.25) 

or 

Sa.  <  <  Sa,  for  5  <  0  (3.26) 

then  phase  transformation  conditions  have  to  be  satisfied  to  find  the  solution  at  time  f  „+i . 
After  checking  which  transformation  is  possibly  active,  the  parameter  A*  is  obtained  by 
enforcing  the  satisfaction  of  the  corresponding  phase  transformation  model.  The  elastic 
right  Cauchy-Green  deformation  tensor  can  now  be  written  using  equation  (3.8). 


=  e*p(-A*el,j9„’;,)Cj;,exp(-A-e5,iV„’;,) 

(3.27) 

where 

JV  - 

l|S»+.ll 

(3.28) 

However,  if  Nn-^i 

commutes  with  C*^j,  taking  logarithms  in  (3.27)  we  get 

symiV„+, 

(3.29) 

which  is  identical  to  the  smedl-strEun  kinematic  relation  that  cam  be  derived  from  (2.39). 


It  is  important  to  note  that  the  steps  proceeding  and  following  the  small-strain  update 
are  purely  kinematic  in  nature,  and  hence,  material  independent.  Thus,  when  it  applies, 


this  procedure  provides  a  material  independent  prescription  for  extending  small-strain 
updates  into  the  finite  deformation  range  within  the  framework  of  miiltiplicative  plasticity. 

The  main  xmderlying  assumption  in  arriving  at  (3.29)  is  that  the  phase  transformation 
flow  direction  commutes  with  the  elastic  predictor  for  large  deformation  as  described 
in  (3.19).  Because  of  the  symmetry  of  this  immediately  requires  that  Nn+i  be 

hkewise  symmetric,  i.e.,  that  the  inelastic  spin  vanish,  and  that  have  the  same 

principal  directions  as  Correspondingly,  n„+i  must  commute  with  as  can  be 

derived  from  (2.37),  for  the  small  strain  formulation.  Using  (2.32)  in  (2.39)  implies  that 
n„+i  commutes  with  if  and  only  if  n„+i  commutes  with  Cn-n  •  By  substituting  elastic 
relations  <r„+i  =  ff’(e^+i)  in  the  phase  transformation  flow  direction  n„+i(tr),  it  can  be 
thought  of  ais  a  function  of  n(e**)  of  the  elastic  strain  Then  the  required  condition 
is  that  the  function  n(e*^)  have  the  same  principal  directions  as  c**,  for  all  vEdues  of  its 
arguments.  A  general  representation  of  a  tensor-valued  function  satisfying  this  requirement 
is 

n(e‘')  =  Ai(e*')/  -I-  A2(e*')e*'  +  A3(c*')c*'’  (3.30) 

where  the  scalar  functions  {Aj,  Aj,  A3},  of  need  not  be  isotropic  (Cuitino  et  al.  [1992]). 
Thus  the  direction  of  the  phase  transformation  flow  n  regarded  as  a  f\mction  of  (e**),  while 
restricted  to  be  of  the  form  (3.30),  cam  still  possess  some  measure  of  anisotropy  (Cuitino  et 
al.  [1992]).  It  bears  emphasis  that  these  requirements  place  no  restrictions  on  the  elastic 
relation  tr„  which  can  be  specified  arbitrarily.  In  particular,  the  elastic 

response  can  be  anisotropic. 

3.5  Boundary  Value  Problem 

In  order  to  formulate  a  well-defined  problem,  in  addition  to  the  constitutive  relations 
we  need  to  consider  the  momentum  balance  equation  and  suitable  boimdary  and  initial 
conditions. 
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3.5.1  Strong  form  of  the  problem 


Let  Qo  he  an  open  set  in  71^  representing  the  image  of  a  body  B  at  time  to-  The 
study  of  the  motion  and  deformation  of  the  physical  body  B  then  reduces  to  the  study 
of  mappings  ;  fio  C  7^*  — »  11^ .  We  assume  a  constitutive  relation  on  ft©  so  that 

the  first  Piola-Kirchhoff  stress  tensor  P  is  a  function  of  P  of  Let  the  deformation 

mapping  be  prescribed  on  the  Dirichlet  boundary  as  #t|a„no  =  The  space  of 
configurations  is  defined  as 


C  =  =  #„  onSho} 

(3.31) 

As  usual  we  require  that 

(3.32) 

dudo  n  deilo  =  0 

(3.33) 

where  5„17o  corresponds  to  the  portion  of  the  boundary  with  prescribed  essential  boundary 

conditions,  while  where  5<rf2o  corresponds  to  the  portion  of  the  boundary  with  prescribed 

natur2d  boundary  conditions. 

The  formed  statement  of  the  strong  form  of  the  problem  is:  Given  a  set  a 

conditions,  find  the  motion  #(X,t)  :  fioxjO,  !r[— »  such  that 

boundary 

DlV{P)  +  pf  =  pU  in  nox]o,r[ 

(3.34) 

=  $  on  5„f2ox]0,r[ 

(3.35) 

P  AT  =  t  on  a,,nox]0,r[ 

(3.36) 

4>(X,0)  =  4»o(X)  allxGfio 

(3.37) 

V(X,0)  =  Vo(X)  allxef2o 

(3.38) 

where  P  is  the  first  Piola  Kirchhoff  stress  tensor,  /  is  the  body  force  per  unit 

mass,  p  is 

the  mass  density,  X  is  the  unit  normal  to  the  boundary,  t  is  the  kirchhoff  stress  vector. 


and  #o(-y)  and  Vb(JC)  are  the  prescribed  initial  displacement  and  velocity  conditions, 
respectively. 

3.5.2  Weak  form  of  the  problem 

The  space  of  admissible  variations  is  defined  as 

V  =  |i7  :  flo  I  *7  =  0  on  (3.39) 

Let  Tf)  =  C  xV  —^Hhe  defined  as 

W»7)  =  /  pA-ijdVo+  f  P.VovdVo-  I  f  vdVo-  f  t  •  d5o  =  0  (3.40) 

Jctfi  •/fio  •/fto 

where  A  =  11.  Since  the  treatment  of  the  transient  dynamic  problem  plays  no  role  in  the 
present  discussion,  we  shall  ignore  inertia  effects  and  confine  our  attention  to  the  static 
case. 

Consider  a  solid  initially  occupying  a  reference  configuration  flo,  and  a  process  of 
incrementally  loading  whereby  the  deformation  mapping  over  flo  changes  from  at 
time  f„,  to  +  ti,  at  time  <„+i  =  <«  +  At.  We  enforce  equilibrium  at  time  <„+i 

weakly  by  recourse  to  the  principal  of  virtual  work.  The  formal  statement  of  the  weak 
form  of  the  boundary  value  problem  is:  Find  #  €  C  such  that  ri)  =  0  for  all  €  V. 

/  Pn+l  :  Vo»7<fVo  —  f  /n+l  •  »ldVo  —  /  ^n+l  ‘  V  =  0  (3-41) 

Jtlo  JQq 

where 

P„+i:  denotes  the  first  Piola-Kirchhoff  stress  field  at  time  <„+i. 

fn+i  ■  is  the  body  force, 

tn-f-i:  are  the  boundary  tractions, 

rj:  is  an  admissible  virtual  displacement  field, 

Vq:  denotes  the  material  gradient. 
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Assuming  for  now  that  we  have  determined  a  rule  to  update  the  stress  field  of  the 
general  form 

Pn+I  =  state  at  t„,  At)  (3-42) 

where  the  deformation  gradient 

Pn+l  =  ^O^n+1  (3.43) 

is  assumed  given. 

The  “consistent  tangents”  K„+i  follow  by  linearization  of  the  stress  updates  (3.42)  as 

dP 

At)  (3.44) 

Orn-i-l 

The  corresponding  “spatial  consistent  tangents”  fc„+i  are  obtained  by  push  forward  to  the 
current  configuration 

kijki  =  KijtiFjjFii  (3.45) 

It  is  important  to  note  that  the  use  of  consistent  tangent  matrices  can  be  computa¬ 
tionally  advantageous  because  of  their  symmetry,  and  that  they  also  have  superior  local 
approximation  properties.  For  small  strains  the  computation  of  the  consistent  tangents 
reduces  to  a  straigntforward  exercise  for  many  commonly  used  models.  However  in  the 

finite  deformation  range  the  calculations  are  considerably  more  cumbersome.  So  in  the 

remaining  of  this  chapter  we  extend  small-strain  updates  and  consistent  tangents  to  the 
finite  deformation  range  by  operating  strictly  at  the  level  of  the  kinematics. 

3.6  Stress  update  strategy 

Within  the  context  of  finite  element  analysis  the  solution  of  problem  (3.40)  is  accom¬ 
plished  by  an  iterative  scheme  such  as  the  Newton’s  method.  Typically  one  solves  a 
sequence  of  linearized  problems  defined  as 

[tr(Vnr  Vu)  +  Vr, :  :  Vu;;i,)]  dV, 

=  -e(<+i.>))  (3.46) 
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until  the  residual  Gi^n+\'>V)  vanishes  to  within  a  prescribed  tolerance. 

The  convergence  rate  of  iterative  scheme  is  essentially  governed  by  the  choice  of  tangent 
moduli  which  depends  in  turn  on  the  iteration  scheme  adopted.  In  a  typical  iteration 
i  +  1  within  a  time  step  (<nitn+il)  the  variables  may  be  obtained  from  either 

(i)  the  values  corresponding  to  the  previous  non-converged  iteration,  or 

(ii)  the  converged  values  (•)„  from  the  previous  time  step. 

Both  procedures  define  algorithms  which  are  consistent  with  the  field  equations.  However, 
scheme  (i)  introduces  a  ‘history  dependence’  of  the  converged  values  on  intermediate  non- 
converged  iterates.  This  may  pose  diflBculties  due  to  the  strong  path  dependence  of  shape 
memory  models.  Spurious  unloadings  at  some  Gauss  points  may  also  occur  as  a  result 
of  this  procedure.  By  contrast,  history  dependence  on  intermediate  non-converged  values 
is  eliminated  with  the  use  of  scheme  (ii),  and  fictitious  numerical  imloading  is  therefore 
prevented.  For  a  detailed  account  of  these  issues  see  Simo  et  al.  [1985]. 

3.7  Derivation  of  the  Consistent  Tangent 

A  lengthy  but  straightforward  calculation  of  the  finite  deformation  spatied  tangents, 
Cuitino  Ortiz  [1992],  gives  the  following  expression. 

^ijkl  ~  ^ijkl  '^nj^inkl  '^ni'^jnkl  *1'  ^ik'^ij 

where  r  is  the  KirchhofF  stress  tensor.  The  various  terms  in  (3.47)  are  defined  as  below. 

<^ijki  =  Fii  Dijmn  Lmnkl{Fii  -I-  FhFij^)  (3.48) 

=  Ft,  Fit  ^  F,d 

=  Fii  Fgl  EjbrsHrsmn  Lmnkl  {Fil  F^f^-  +  Ff^-)  (3.49) 

where  D  represents  the  samll-strain  tangent  defined  as 


which  is  the  only  material  dependent  quantity  in  eq  (3.47).  F‘  may  be  thought  of  as  an 

elastic  predictor  for  F*,  and 

(3.51) 

are  the  tangent  elastic  compliances. 

In  component  form 

f:;  =  FAiF^-'uj 

(3.52) 

j.  d{\og\/^)MN 

^ M N Lt 

E  - 

a(A<eJ,JV)«j 

(3.53) 

a[cxp(A*  e*i^N)]ij 

„  d{XUiN)Rs 

^RSMN  =  a  e* 

(3.64) 

=  {^RM^SN  +  ^RN  ^Sm)  -  C%sAB  ^ABMN 

(3.55) 

Crsab  -  ha)  3^  ^  ^^((ba  (rs) 

(3.56) 

where  A  and  n  in  the  last  equation  are  the  Lame’s  constants. 

3.7.1  Derivation  of 

,  a(iog\/^) 

^  ac* 

(3.57) 

By  the  spectral  theorem  for  symmetric  positive  definite  tensors  we  have 

3 

C'  =  J]  0  iV^“>  (o  =  1, 2, 3) 

a=I 

(3.58) 

The  invarients  of  C*  are  related  to  the  coefficients  in  the  characteristic  polynomial  p(A^) 

of  C*  as  follows, 

p{Xl):=-Xl  +  hXl-l2Xl  +  h  =  0  (a  =  l,2,3) 

(3.59) 
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where  >  0  are  the  roots  of  the  chiiracteristic  polynomial.  Therefore 

3 

a=l 

3 

=>  log\/^  =  ^  logA„iV(“)  ®  Ar(“)  (3.60) 

a=l 

Note:  If  is  a  proper  vector  of  C*,  it  is  also  a  proper  vector  of  log\/ C®.  The  principal 
invarients  of  C*  6  are  denoted  by 

h  ;  =  tr(C*];  /2  ^  [/?  -  tr(C'')]  ;  1,  :=  dei[C‘]  =  (3.6I) 

where 

5^  =  |C*  G  Afi  |C'  =  C*’’}  (3.62) 

Following  standard  usage,  we  denote  by  'R?)  the  vector  space  of  linear  transformation 

in  TV,M.^  is  defined  as 

M\  =  {F€r(7^^:7^^)|det(F]  >0} 

By  the  chain  rule 

51og\/^  51og\/^  dl\  d\ogy/^  5/2  91og\/^  dl^ 

dc  ^  a/I  dc  a/2  ^ '''  ais  ^ 

In  order  to  evaluate  the  above  expression  we  need  the  following  results. 

aiog\/^  aiog\/i^  aAi  aiog\/^  aA2  aiogv/c*  aAs 

dia  ~  ~~dXi  dl^  dX2  5^  ~dX^  dh 

3. 7. 1.1  Proposition:  Assume  that  Aj  ^  A2  ^  A3.  By  (3.59)  Aa  is  a  function  of  the  invariants 
of  C*,  denoted  by 

A.  =  A.U„/2,/3)  (»  =  1,2,3) 

Then  the  following  relations  hold: 

dXa  _  XI  dXa  -A„  dX,  A;2 

a/,  2Da'  dh  2D  a'  a/3  Da 
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(3.66) 


(3.67) 


(3.63) 


(3.64) 


(3.65) 


where  Da  =  (Aj  —  Aj)(A^  —  A^).  See  proof  in  Simo  et  al.  [1991],  p.  227. 
Consequently,  we  can  derive  the  following  relations  as  follows 


^logv^  =  A^iogA„7V<“)®iV<“) 


dh 


dh 


a=l 


Similarly,  for  the  second  term  we  get 


“) 


and 


d 

dh 


log  =  I J2  (^)  ® 


Using  the  chain  rule  for  the  second  term  in  equ  (3.64) 

dia  _  dia  5Ai  dia  d\2  dia  8X3 

dC‘  dXi  ac*  dXi  dC‘  6X3  dC‘ 


where 


^^°=2vKt  (;9=  1.2,3) 


aXg 


(3.68) 


(3.69) 


(3.70) 


(3.71) 


(3.72) 


3. 7. 1.2  Proposition:  Let  the  spectral  decomposition  of  C'  be  defined  by  eq.  (3.58)  So  Aa 
is  a  fimction  of  C*  and  we  write  Aa  =  Aa(C*).  Upon  defining 


we  have  the  basic  relations 


(3.73) 


(1)  ^  (<■  =  1.2.3)  iffA, /A2#A. 

(2)  ^  =  jA.lC--’  -  M<’»]  a  A,  =  Aj  5^  A3 

^  =  jAsMC) 

ac*  2 

(3)  =  JaiC‘'‘ 

ac*  2 


iff  Aj  =  A2  =  A3 


Consequently,  using  these  relations  in  (3.71),  we  get 

3 


ac 


^  =  E21‘  (5A.M<‘>) 

=  (At-2iV(‘)(8)Ar(‘)) 

t=i 

=  Z!  ® 


6=1 


Similarly,  for  the  second  term  we  get 


a/2 


ac 


7  =  E 


6s  1 
3 


a/2  aAi 

aAfc  ‘  ac* 


=  E2A,{/,-A5)-Qa,MI‘>) 
=  Y,{h-  Ai)  iV<*)  ® 


6=1 


and 


a/3 


a/3  dXb 


ac*  ~  H  dXi, '  ac* 


6=1 

3 


(3.74) 

(3.75) 

(3.76) 

(3.77) 


(3.78) 


(3.79) 


(3.80) 


B=l 


Remark:  For  a  general  framework  see  Marsden  &  Hughes  [1983],  pp.  220-222. 
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Substituting  the  corresponding  expressions  in  (3.64) 

a(iog%/o) 


ac« 


+ 


+ 


am 


J  (8)  ■ 

•  (7i  -A^)Ar(‘)0N(‘> 

\fc=i 

)  0  • 

where  D.  =  (Aj  -  Aj)(Aj  -  A?) 

The  closed  form  expression  for  is 


3. 7. 1.3  First  Order  Approximation 

A  first  order  approximation  for  v/^  can  be  expressed  as 

Consequently,  a  first  order  approximation  for  logV C'  c£in  be  written  as 
log\/3^  5  log  Ql/  +  C'j)  =  1(C«  -  /)  -  1(C'  -  Jf  + 


Therefore 


5(log\/^)  MN 


dC 


KL 


dCK 


/)mA.  -  -  I)\,N 


} 


3  1  - 

=  ^  fiMK  ^NL  -  g  +  Sf^i 


(3.81) 


(3.82) 


(3.83) 


(3.84) 


(3.85) 
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3.7.2  Derivation  of  Hrsmn 


^  a(exp(VeiJV(5)))  . 

(X-eiNiS))  " 


(3.86) 


Let  (A*e^  A/’(S))  =  K.  Employing  the  spectral  theorem  for  symmetric  positive  '^“finite 
tensors,  (3.58),  we  can  write 


exp(Ji^’)  =  ^exp(Aa)N^°^  0  N^‘ 


(3.87) 


Aa’s  in  (3.88)  should  not  be  confused  with  previously  defined  Aa’s.  The  principad  invarients 
of  Jf  G  5^  are  denoted  by 


/,  :  =  trfAT];  I,:=\  [/?  -  tr(A:2)]  ;  1,  :=  det[K]  = 


(3.88) 


where  5|  =  {K  €  Ml\K  =  K'^}  We  denote  by  C{TV,7l^)  the  vector  space  of  linear 
transformation  in  71^,  is  defined  in  eq.  (3.63).  Therefore 


5exp(K)  _  5exp(J^)  dli  dexp{K)  dh  dexp{K)  dh 

dK  ~  dh  dK  ~~dh  ^  dh  ^ 

Once  again,  in  order  to  evaluate  the  above  expression  we  need  the  following  results. 


(3.89) 


dexp(K)  _  d€xp(K)  dXi  ,  dexp(K)  6X2  ,  dexp{K)  dXj 

dia  dXi  dl^  dX2  dh  dX3  dh  ^  ’ 


Assuming  that  Aj  ^  A2  ^  A3.  By  (3.59)  A^  is  a  function  of  the  invariants  of  K,  denoted 
by 

Aa  =  Aa(/i,/2,/3)  (a  =  1,2,3)  (3.91) 

Then  the  following  relations  hold: 

-  M!.  ' 

dl,  ■  D, '  dh  D.  '  dh  D. 
where  Da  is  defined  earlier.  For  proof  see  Simo  et  al.  [1991],  p.  227. 
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Consequently, 


Similarly, 


and 


^exp(Jif)  =  A  j2exp(A;)JVl“)®ArCl 

^exp(J«:)  =  |;exp{A;)(^)Nl">®JV<“> 
Aexp(Jir)  =  ^mp(AJ)  iVl“>  ®  N'”' 


Employing  the  chaiin  rule  for  the  second  term  in  (3.90), 

dia  dia  d\\  dh  d\l  dh  d\l 
dK  ~  dX]  dK  dXl  dK  ^  dX\  dK 


where 


dA^  Q\^ 
dK  ~  W^'dK 

dK 


Substituting  in  (3.90),  we  get 
5(explif) 


(3.92) 


(3.93) 


(3.94) 


(3.95) 


(3.96) 


+  ^J^exp(A2j(z^)  .  ^^(/,  -A2)iV<'^®N<‘>^ 

+  ^^exp(A2)^^^fV<“)®  •  ^^(/3Aj-')ArW®iV<‘>^  (3.97) 

where  Da  =  (A„  -  Aj)(Ao  -  A^)  Using  (97)  in  (86)  we  get  the  desired  result. 


3. 7.2.1  First  Order  Approximation 


In  component  fonn,  we  can  write  exp(Kij)  as  follows 

exp( X.j)  a  fij  +  ATy  +  +  jfK.iKuKij  +  ■■■  (3  98) 

Retaining  up  to  quadratic  terms  in  (3.99)  we  get 

+  K,kSji)  (3.99) 

3.7.3  Finite  Deformation  Update 

Once  all  terms  in  eq.  (3.48)  and  (3.49)  are  calculated,  the  full  finite  deformation  update 
comprises  three  steps: 

1.  Preprocessor  step.  Compute  the  predictor  logarithmic  elastic  strains  jlogC*^j  from 
the  given  updated  deformation  gradients  as  in  (3.19).  Set  equal  to  jlogC'^j . 
Identify  <r„  with  S„. 

2.  Small  strain  update.  Effect  a  small-strain  update  driven  by  c**^i  with  initial  condi¬ 
tions  a„,  to  compute  <r„+i  and  A*. 

3.  Postprocessor  step.  Identify  .Sn-t-i  with  and  compute  exponentiaJ 

mapping  (3.8). 

It  is  important  to  note  that  the  steps  preceding  and  following  the  small-strain  update 
are  purely  kinematic  in  nature,  and  hence,  material  independent. 


3.8  Numerical  Results 


In  order  to  test  the  integration  algorithm  in  the  finite  deformation  range,  we  performed 
a  series  of  imiaxial  simulations  \mder  displacement  control.  The  material  properties  are 

E  =  7500kPa,  »/ =  0.3  (3.100) 

while  the  transformation  ranges  are: 

om,  =  50,  a  Ms  =  70,  a  A,  =  75,  a  a,  =  90  (3.101) 

with 

Cm  =  1,Cx  =  1,  cl  =  0.6  (3.102) 


The  tests  are  organized  as  follows: 

•  Test  1;  we  run  a  loading- unloading  test  with  complete  A  M  transformation  and 
complete  M  —*  A  transformation  (see  Figures  3.8.1  -  3.8.5). 

•  Test  2:  we  run  a  loading- unloading  test  with  incomplete  A  —*  M  transformation  and 
complete  M  —*  A  transformation  (see  Figures  3.8.6  -  3.8.10). 

•  Test  3:  we  run  a  loading-unloading  test  with  complete  A  —*  M  transformation  and 
incomplete  M  —*  A  transformation  (see  Figures  3.8.11  -  3.8.15). 

•  Test  4:  we  run  a  loading- unloading  test  with  incomplete  A  —*  M  transformation  and 
incomplete  M  —*  A  transformation  (see  Figures  3.8.16  -  3.8.20). 


3.8.1  Test  1:  Complete  A  M  and  M  A  transformations 

Figure  3.8.1  presents  the  stress-strain  plot  for  complete  phase  transformations  in  either 
case.  The  applied  strain  as  a  function  of  the  time  parameter  is  presented  next.  Veiriable 
loading  increments  can  be  applied  under  the  present  technique.  The  specimen  is  loaded  to 
a  certain  strain  and  imloaded  to  zero  stress  condition  in  each  loading  cycle.  Figure  3.8.3 
presents  the  stress  variation  at  a  Gauss  point  as  a  function  of  t.  Figures  3, 8.4-5  show  the 
variation  in  the  internal  phase  transformation  parameter  ^  as  a  fimction  of  strain  e  at  the 
Gauss  point  and  of  time  (<),  respectively. 


Pseudo-elastic  Behavior 

Complete  M  — >  A,  Complete  A  — >  M 


Fig.  3.8.1  Stress  Vs  strain  plot.  Complete  M  A  Eind  complete  A  —*  M  transformations. 


Fig.  3.8.2  Strain  Vs  time  plot.  Complete  M  —*  A  and  complete  A—*M  transformations. 


Fig.  3.8.3  Stress  Vs  time  plot.  Complete  M  A  and  complete  A  M  transformations. 
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Fig.  3.8.4  Internal  parameter  Vs  strain.  Complete  M  —*  A  and  complete  A  —*  M  transformations. 


Fig.  3.8.5  Internal  parameter  Vs  time.  Complete  M  A  and  complete  A  M  transformations. 
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3.8.2  Test  2:  Incomplete  A-*  M  and  complete  M  A  transformations 

Figure  3.8.6  presents  the  stress-strain  plot  for  incomplete  A  —*  M  and  complete  M  —*  A 
transformation.  The  applied  strain  as  a  function  of  the  time  parameter  is  presented  next. 
The  strain  is  increased  gradually  and  the  specimen  is  unloaded  to  zero  stress  condition  in 
each  cycle.  Figure  3.8.8  present.s  the  stress  variation  at  a  Gauss  point  as  a  function  of  t. 
Figures  3.8.9-10  show  the  variation  in  the  internal  phase  transformation  parameter  ^  as  a 
function  of  strain  e  at  the  Gauss  point  and  of  time  (<),  respectively. 


Pseudo-elastic  Behavior 

Complete  M  — >  K  Incomplete  A  — >  M 


Fig.  3.8.6  Stress  Vs  strain  plot.  Complete  M  A  and  incomplete  A  M  transformations. 


Fig.  3.8.8  Stress  Vs  time  plot.  Complete  M  —*  A  and  incomplete  A  M  transformations. 


Fig.  3.8.9  Internal  parameter  Vs  streiin.  Complete  M  A  and  incomplete  A  —*  M  transformations. 


Fig.  3.8.10  Internal  parameter  Vs  time.  Complete  M  -*  A  and  incomplete  A  —*  M  transformations. 
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3.8.3  Test  3:  Complete  A —*  M  and  incomplete  M  A  transformations 

Figure  3.8.11  presents  the  stress-strain  plot  for  complete  A  —*  M  and  incomplete 
M  —*  A  transformation.  The  applied  strain  as  a  function  of  the  time  parameter  is  presented 
next.  The  specimen  is  imloaded  to  zero  stress  condition  at  the  end  of  the  final  loading 
cycle.  Figure  3.8.13  presents  the  stress  variation  at  a  Gauss  point  as  a  function  of  t. 
Figures  3.8.14-15  show  the  variation  in  the  internal  phase  transformation  parameter  ^  as 
a  function  of  strain  e  at  the  Gauss  point  and  of  time  (<),  respectively. 


Pseudo-elastic  Behavior 

Incomplete  M  — >  A,  Complete  A  — >  M 


Fig.  3.8.11  Stress  Vs  strain  plot.  Incomplete  M  —*  A  and  complete  A  M  transformations. 
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Fig.  3.8  .12  Strain  Vs  time  plot.  Incomplete  M  -*  A  and  complete  A  M  transformations. 
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Internal  Parameter  (xi) 


Fig.  3.8  .14  Internal  paurameter  Vs  strain.  Incomplete  M  -*  A  and  complete  A—*M  transformations. 


Fig.  3.8.15  Internal  parameter  Vs  time.  Incomplete  M  A  and  complete  A  M  transformations. 


3.8.4  Test  4:  Incomplete  A  —*  M  and  incomplete  M  —*  A  transformations 

Figure  3.8.16  presents  the  stress-strain  plot  for  incomplete  A  —*  M  and  incomplete 
M  —*A  transformations.  A  random  strain  loading  is  applied  as  a  function  of  the  time 
parameter  jmd  the  specimen  shows  a  strong  path  dependence,  depicted  in  Fig.  3.8.16.  Figure 
3.8.18  presents  the  stress  variation  at  a  Gauss  point  as  a  function  of  i.  Figures  3.8.14-15 
show  the  variation  in  the  internal  phase  transformation  parameter  ^  as  a  function  of  strain 
e  at  the  Gauss  point  and  of  time  (<),  respectively. 


Pseudo-elastic  Behavior 

Incomplete  M  — >  A,  Incomplete  A  — >  M 


Fig.  3.8.16  Stress  Vs  strain  plot.  Incomplete  M  —*  A  and  incomplete  A  —*  M  transformations. 
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Fig.  3.8.17  Strain  Vs  time  plot.  Incomplete  M  A  and  incomplete  A  —*  M  transformations. 
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Fig.  3.8  .19  Internal  parameter  Vs  strain.  Incomplete  M  A  and  incomplete  A  M  transformations. 


Fig.  3.8  .20  Internal  parameter  Vs  time.  Incomplete  M  —*  A  eind  incomplete  A  —*  M  transformations. 


CHAPTER  4 


A  FINITE  ELEMENT  THEORY  FOR  COMPOSITE  LAMINATES 
SPECIALIZED  TO  FLAT  GEOMETRIES 


4.1  Introduction 

In  the  last  two  decades,  composites  have  foimd  increasing  application  in  many  engineer¬ 
ing  structures.  Recent  advances  in  the  technologies  of  manufacturing  and  materials  have 
enhanced  the  current  application  of  composite  materials  from  being  used  as  secondary 
structural  elements  to  becoming  primary  load-carrying  structural  components.  Due  to 
the  inherent  inhomogeniety  suid  anisotropy  of  the  materials,  analysis  of  these  composite 
structures  imposes  new  challenges  to  engineers  (Noor  &  Burton  [1989,1990]). 

Plate  and  shell  structures  made  of  laminated  composite  materials  are  often  modeled 
as  an  equivalent  single  layer  using  classical  laminated  theory  (C.L.T)  (Christensen  [1979], 
Jones  [1975],  Pagano  [1989],  Pipes  &  Pagano  [1970])  in  which  the  thickness  stress  com¬ 
ponents  are  ignored.  The  classical  laminate  theory  is  a  direct  extension  of  classical  plate 
theory  in  which  the  well  known  KirchhofF-Love  kinematic  hypothesis  is  assumed  enforced. 
This  theory  is  awlequate  when  the  thickness  (to  side  or  radius  ratio)  is  smeJl.  However, 
laminated  plates  and  shells  made  of  advanced  fUamentry  composite  materials  are  succep- 
tible  to  thickness  effects  because  their  effective  transverse  moduli  are  significantly  smaller 
than  the  effective  elastic  modulus  along  the  fiber  direction.  Furthermore,  the  classical 
theory  of  plates  which  assvunes  that  the  normals  to  the  midplane  before  deformation  re¬ 
main  otraight  and  normal  to  the  plane  eifter  deformation,  imderpredicts  deflections  and 
overpredicts  natural  frequencies  and  buckling  loads.  These  discrepencies  are  due  to  the 
neglect  of  transverse  shear  strains.  The  errors  in  deflection,  stresses,  natiiral  frequencies, 
and  buckling  loads  are  even  higher  for  plates  made  of  advanced  composites.  The  range 


of  applicability  of  the  C.L.T.  solution  has  been  well  established  for  laminated  flat  plates 
(Noor  [1989],  Pagano  [1989],  Reddy  [1993]).  These  analyses  have  indicated  that  a  theory 
which  accounts  for  the  transverse  shear  deformation  effects  would  be  adequate  to  predict 
the  gross  behavior  of  the  laminate. 

In  order  to  overcome  these  deficiencies  in  C.L.T. ,  refined  laminate  theories  have  been 
proposed  (Byun  &  Kapania  [1992],  Noor  k  Burton  [1989],  Reddy  [1984,  1993],  Robbins 
et  al.  [1991],  Whitney  [1970]).  These  are  single  layer  theories  in  which  the  transverse 
shear  stresses  are  taken  into  account.  They  provide  improved  global  response  estimates 
for  deflections,  vibration  frequencies  and  buckling  loads  of  moderately  thick  composites 
when  compared  to  the  classical  laminate  theory.  A  Mindlin  type  first-order  transverse 
shear  deformation  theory  (S.D.T.)  was  first  developed  by  Whitney  and  Pagano  [1970]  for 
multilayered  anisotropic  plates,  and  by  Dong  and  Tso  [1972]  for  multilayered  anisotropic 
shells.  Both  approaches  (C.L.T.  and  S.D.T.)  consider  all  layers  as  one  eqmvalent  single 
anisotropic  layer;  thus  they  are  inadequate  to  model  the  warping  of  cross-sections,  that 
is,  the  distortion  of  the  deformed  normal  due  to  transverse  shear  stresses.  Furthermore, 
the  assumption  of  nondeformable  normeil  results  in  incompatible  shearing  stresses  between 
adjacent  layers.  The  later  approach  also  requires  the  introduction  of  an  arbitrary  shear 
correction  factor  which  is  dependent  on  the  lamination  parametrs  for  obtaining  accurate 
results. 

The  exact  analyses  performed  by  Pagano  [1970,1973,1989]  on  the  composite  flat  plates 
have  indicated  that  the  distortion  of  the  deformed  normjil  is  dependent  not  only  on  the 
laminate  thickness,  but  also  on  the  orientation  and  the  degree  of  orthotropy  of  the  indi¬ 
vidual  layers.  Therefore  the  hypothesis  of  nondeformable  normals,  while  acceptable  for 
isotropic  plates  2ind  shells  is  often  quite  unacceptable  for  multilayered  anisotropic  plates 
and  shells  with  very  large  ratio  of  Young’s  modulus  to  shear  modulus,  even  if  they  are 
relatively  thin  (Chang  et  al.  [1990],  Reddy  [1989,  1993]).  Thus  a  transverse  shear  defor- 
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mation  theory  which  also  accounU  for  the  distortion  of  the  deformed  normal  is  required  for 
accurate  prediction  of  the  elastic  linear  behavior  (deflections,  thickness  distribution  of  the 
in-pliine  displacements,  natural  frequencies,  etc.)  of  multilayered  anisotropic  plates  and 
shells. 

In  view  of  these  issues  a  variationally  soimd  theory  that  accounts  for  the  3-D  effects, 
allows  thickness  variation,  and  permits  the  warping  of  the  deformed  normal  is  required  for 
refined  and  accurate  analysis  of  thick  and  thin  composites.  Embedding  of  shape  memory 
alloys  and  study  of  its  interaction  with  composite  laminates  only  in  the  context  of  such  a 
theory  will  be  feasible. 

4.2  New  Ideas  Proposed  in  the  Present  Theory 

This  section  presents  the  noteworthy  attributes  of  the  proposed  theory  that  directly 
address  some  of  the  technical  drawbacks  present  in  most  of  the  finite  element  theories  that 
have  been  proposed  for  composite  analysis. 

(i)  The  displacement  field  is  continuous  through  the  thickness  of  the  composite  while  the 
rotation  field  is  layer-wise  continuous  (in  2-D)  and  can  be  discontinuous  across  the 
finite  element  layers  through  the  thickness  direction.  This  displacement  field  fulfills  a 
priori  the  geometric  continuity  conditions  between  contiguous  layers. 

(ii)  The  assumed  displacement  field  is  capable  of  modeling  the  distortion  of  the  deformed 
normal,  without  increasing  the  order  of  the  peirtial  differential  equations  with  respect 
to  the  first -order  transverse  shear  deformation  theory. 

(iii)  The  assumed  displacement  field  has  a  3-D  feature,  thereby  modeling  accurately  the 
interlaminar  conditions  and  predicting  the  3-D  edge  effects. 

(iv)  Like  the  shear  deformable  theory  the  proposed  composite  shell  theory  provides  flexi¬ 
bility  in  the  spacification  of  the  boundary  conditions,  Hughes  [1987]. 

(v)  In  this  theory,  at  most,  only  first  derivatives  of  displacement  Eind  rotation  fields  appear 
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in  the  variational  equations.  The  practical  consequence  of  this  fact  is  that  only 
continuity  of  finite  element  functions  is  required  which  is  readily  satisfied  by  the  family 
of  Lagrange  elements,  Hughes  [1987]. 

(vi)  It  is  feasible  to  employ  this  formulation  for  constructing  plate  and  sheU  finite  elements 
via  the  finite  element  displacement  method,  Hughes  [1987]. 

4.3  Assumptions  of  the  Layer-wise  Shear  Deformable  Shell  Theory 

1.  The  domain  fl  is  of  the  following  special  form 

x,yf^  €  C 

where  is  the  area  of  the  reference  surface  for  layer  ‘/’  and  T  is  the  total  thickness 
of  the  composite  shell. 

2.  The  displacement  field  is  assumed  to  take  the  following  form 

Vi'\x,y,z)  =  +  (4.2) 

where  Ua*(x,y,  z)  are  the  translations  and  da\i,y)  are  the  out  of  plane  director  rota* 
tions  for  the  /**  layer  reference  surface,  and  C  6]0, 1[  is  a  parameter  that  establishes 
the  position  of  a  point  from  the  reference  surface. 

3.  The  displacement  field  in  the  thickness  direction  is  assumed  to  be  a  function  of  z,  i.e., 

^3*^  =  «3(a^,y,2)  (4.3) 

and  by  virtue  of  this  assumption  the  transverse  displacement,  U3,  does  vary  through 
the  thickness,  thereby  producing  through  the  thickness  strains  which  result  in  thickness 
variation  in  the  shell. 

4.  As  a  consequence  of  the  above  relajcation 

(4.4) 


(4.1) 


z)  e 


o  ’ o 


<^33  /  0 
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i.e.,  we  do  not  invoke  the  plane  stress  hypothesis. 


5.  We  know  from  the  elasticity  theory  that  the  displacements  and  stresses  at  the  interface 
between  and  (/  +  1)*‘  bounded  layers  must  satisfy  the  following  contact  conditions. 

£,(1)  =  (4,5) 

=  Ui'*"  (4.6) 


In  addition,  the  following  stress  conditions  need  to  be  satisfied 


Jf)  _  J'+i) 

^a3  “  ^a3 


Ji+i) 

^33 


(4.7) 

(4.8) 


where  denotes  the  interlaminar  stress  for  layer  /. 

A  consequence  of  the  second  assumption  is  that  each  finite  element  layer  is  associated 
with  non-normal  cross-sectional  rotations  in  accordance  with  the  Mindlin  kinematic  as¬ 
sumption.  Another  consequence  of  the  second  assumption  is  that  it  results  in  independent 
sheeir  deformation  of  the  director  in  each  layer  and  allows  the  warping  of  the  composite 
cross-section.  It  also  results  in  discontinuous  strain  fields  across  the  different  material  sets, 
thereby  creating  the  provision  of  stress  continuity  across  the  materied  interfaces.  Conse¬ 
quently,  equations  (4.7)  and  (4.8)  are  inherently  satisfied. 


4.4  Kinematic  Description  of  Multi-layered  Shells 

Figure  4.1  shows  a  composite  laminate  with  ‘m’  number  of  material  layers.  (In  Fig.  4.1, 
m  =  3).  Consequently,  the  total  thickness  T  is  divided  into  ‘m’  finite  elements  through 
the  thickness.  Each  layer  of  finite  elements  is  associated  with  a  reference  surface  which, 
for  ease  of  discussion,  is  assumed  to  be  coincident  with  the  lower  surface  of  that  layer. 
It  is  important  to  note  that  the  director  rotation  (i.e.,  ^q)  and  the  slope  (i.e.,  us  o)  are 
not  necessarily  the  same  and  thus  transverse  shear  strains  are  accommodated.  This  is  to 
be  contrasted  with  the  classical  lamination  theories  (C.L.T.)  in  which  $q  =  U3  ^  and 
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therefore  the  transverse  shear  strains  are  zero.  Within  each  layer  (of  finite  elements), 
the  director  rotates  by  0a  \  generating  shesir  strain  Consequently,  in  the  deformed 
configuration,  node  2  (see  Fig.  4.1)  moves  to  location  2'.  Now,  the  director  in  the  second 
layer  rotates  by  an  angle  0a generating  shear  strain  7^^*^  This  kinematics  repeats  in 
the  successive  layers  said  due  to  the  continuity  of  the  displacement  field  in  the  thickness 
direction  we  obtain  the  new  locations  of  finite  element  nodal  points  as  1  ,  2  ,  3  and  4  . 
This  new  location  of  points  produces  a  higher  order  variation  of  displsu;ement  field  through 
the  thickness  direction,  and  thus  precludes  the  need  for  introducing  arbitrary  polynomial 
expressions  to  model  the  higher  order  variation  of  displacement  field  through  the  thickness. 


Figure  4.1  Shell  kinematics.  Vso’iation  of  transverse  shear  strains  through  the  thickness. 


4.4.1  Kinematic*  in  the  context  of  finite  element  method 


The  displacement  of  the  shell  is  defined  by  the  following  relations 


•?,C)  =  »">({.'))  + 

(4.9) 

asl 

(4.10) 

asl 

(4.11) 

(no  sum) 

(4.12) 

where 

:  is  vhe  displacement  of  a  generic  point  in  the  shell  layer  /, 

;  is  the  displacement  of  a  point  on  the  reference  surface  of  the  shell  layer  I, 
l/(')  :  is  the  ‘director  displacement’  for  the  shell  layer  /,  and 

ZU) 

:  is  the  thickness  function. 

The  director  nodaJ  displacement  vector  is  constructed  so  that  the  director  can  rotate 
and  stretch,  shown  by  an  additive  decomposition  as  follows 

+  Ui'^  (4.13) 


The  tlirough  thickness  stretch  component  can  be  expressed  as 


=  («,3 


(0/ 

a3 


(4.14) 


where  Uas  and  M(a+n,„)3  ^re  the  translations  in  the  thickness  direction  of  the  lower  and 
upper  surfaces  belonging  to  layer  The  out  of  plane  rotation  component  of  the  director 
is  constructed  such  that  it  may  rotate,  viz.. 


(0 

a2 


gU)  (/)/ 
•'ol  ®o2 
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(4.15) 


The  quantities  6^^  and  ffj/2  represent  the  rotations  of  the  director  about  the  basis 
vectors  and  ,  for  shell  layer  /,  respectively. 

4.5  Geometric  Description  of  Multi-layered  Shells 

Figure  4.2  shows  a  typical  configuration  of  a  doubly  curved  composite  shell.  It  can 
be  made  of  numerous  plies  with  variable  material  properties,  ply  orientations  and  ply 
thicknesses.  In  the  proposed  modeling  of  the  mechanics  of  shell,  no  such  eissumption 
hsis  been  made  which  would  limit  the  number  of  individual  plies,  ply  thicknesses,  their 
orientation  or  their  stacking  sequence. 


Figure  4.2  Micro  and  macro  structure  of  the  composite. 

Let  us  concentrate  our  attention  to  the  case  in  which  the  composite  is  composed  of  four 
laminates  (see  Fig.  4.3),  each  being  modeled  via  a  finite  element  layer.  Figure  4.3  shows  a 
schematic  diagram  of  the  geometry  of  an  individual  layer 
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Figure  4.3  Geometric  description  of  an  individuEil  layer. 


4.5.1  Geometric  description  in  the  context  of  the  finite  element  method 

The  geometry  of  a  typical  Isiminate  of  the  composite  shell  element  is  defined  by  the 
following  relations 


(4.16) 

*"’({.•?,«  =  (4.17) 

asl 

JC*''(£.-),0  =  ^NM.v)Xi‘\0  (4.18) 

Xi'HO  =  (no  sum)  (4.19) 


:  the  position  vector  of  a  generic  point  of  the  shell  for  layer  I, 

:  the  position  vector  of  a  point  in  the  reference  surface  for  layer  I, 

:  a  position  vector  of  a  generic  point  relative  to  which  defines  the  director 

through  the  point  for  layer  /.  (In  computational  shell  literature,  X^^^  is  referred  to  as 
fiber  direction), 
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Xa  :  the  position  vector  of  nodal  point  *a’  in  layer  /, 

Na  :  two-dimensional  shape  function  associated  with  node  ‘a’, 
rten  •  the  number  of  element  nodes  in  the  reference  surface  of  layer  I, 

Xa^  :  a  unit  vector  emanating  from  node  ‘a’  in  the  director  direction, 

Za  •  a  ‘thickness  function’  associated  with  node  ‘a’,  which  is  defined  by  the  location 
of  the  reference  surface. 


The  above  relations  represent  a  smooth  mapping  of  the  biunit  cube  into  the  physical 
shell  domain.  For  ‘C’  fixed,  the  surface  defined  by  (4.16)  is  called  a  lamina  and  for 
fixed,  the  line  described  by  (4.16)  is  called  the  director.  The  dire  ctor  is,  in  general,  not 
perpendicular  to  the  laminae.  Sometimes  the  director  is  referred  to  as  the  ‘pseudonormal’. 


For  a  particular  choice  of  two-dimensional  shape  functions,  eqs.  (4.16)-(4.19)  are  pre¬ 
cisely  defined  upon  specification  of  Xa\  Xa\  and  Za^~  (a  =  1,2,.. .r»en)-  If  i® 

convenient  in  practice  to  take  as  input  the  coordinates  of  the  top  and  bottom  surfaces  of 
the  shell  layer  along  each  nodal  director  {xa^*  and  *9^",  respectively)  and  a  parameter 
(  €  [— 1,-f-l],  which  defines  the  location  of  the  reference  surface.  For  example  if  C  = 
-1,0, -fl  (respectively),  then  the  reference  surface  is  taken  to  be  the  bottom,  middle,  top 
(respectively)  of  the  shell  layer  /.  FVom  these  data  we  may  calculate 

+  (4-20) 


(4.21) 


(C)  =  5(1  +  +  5(1-0  Zy  (4.22) 


where  ||  •  ||  denotes  the  Euclidean  norm  (i.e.,  |j®l]  =  (jj  +  +  ij)'/^). 


4.5.2  Lamina  coordinate  system 

At  each  integration  point  in  the  element  a  Cartesian  reference  frame  is  erected  so 
that  two  axes  are  tangent  to  the  lamina  through  the  point.  The  frame  is  defined  by  its 


orthonormal  basis  vectors  ef,  63 , 63  in  which  63  is  pe  endicular  to  the  lamina.  The  basis 


Figure  4.4  Typical  lamina  coordinate  system  (C  =  constant  surface) 


Construct  unit  tangent  vectors  to  the  and  rf—  coordinate  directions; 


= 


4'>  = 


Consequently 


«3  - 


X  c 


(0 


lie 


'4  ^  e»»  II 

The  vectors  tangent  to  the  lamina  are  selected  so  that  the  angle  between  ef  and  is  the 


(0i 


(4.23) 

(4.24) 

(4.25) 


same  as  the  angle  between  and  Cj  and  so  that  the  ef,  C2 -basis  is  as  ‘close’  as  possible 
to  the  basis.  Thus 


ev) 


4  =  f  (4>  +  e|,") 


(4.26) 

(4.27) 


where 


5(4'^  + 


1 5  (4  ’  + 
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(4.29) 


ea  X  cg^ 

Ca  ^  *0^1 


Figure  4.5  rj  =  constant  surface. 

Furthermore,  we  also  define  the  orthogonal  matrix  to  transform  quantities  from  the 
global  coordinate  system  to  the  lamina  system 


9  =  (9.j1  =  (cf 


(4.30) 


For  the  present  case 


cf  =  Ai  ,  =  D 


(4.31) 


4.6  Constitutive  Relations 


Figure  4.6  Schematic  diagram  of  a  composite  and  its  constituent  laminates. 
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Suppose  the  composite  is  made  of  two  laminates  of  the  same  material  but  oriented 
at  -hff  and  —0  degrees  with  respect  to  the  extension  direction.  Let  C(o)  represent  the  C 
matrix  (constitutive  matrix)  for  the  laminate  with  regard  to  its  mutually  perpendiculEir 


planes  of  elastic  symmetry.  In  general  it  can  be  written  as 
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where  q  represents  the  laminate  in  the  composite. 


(4.32) 


Let  Xj  represent  the  direction  along  the  loading  for  the  composite  element  and  let  X3 
represent  its  thickness  direction.  This  axis  is  assumed  to  be  perpendicular  to  the  plane  of 
elastic  symmetry.  The  C(a)  for  a  laminate  can  be  projected  from  its  mutually  perpendic¬ 


ular  planes  of  elastic  symmetry  onto  the  composite  coordinate  system  (Xi,X2,X3)  about 


the  X3  axis  via  the  following 

transformation  matrix. 
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(4.33) 


where  C  =  cos  0,  S  =  sin  9.  The  transformed  constitutive  matrix  is  obtained  as 


^(q)  -  Qjc,)^(a)Q(Q) 
In  general,  will  have  the  following  form 
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We  evaluate  the  C(„,acro  e/em.)»  i  the  constitutive  relation  for  the  composite  via  the 


expression 


^(macro  eUm.)  —  *^^(1)  "i"  ^)  ^(: 


where  m  =  (1/(2  <  1,  (6  =  ■?(!),  6  =  ■f(l)  +  ■?(2)). 


4.7  Finite  Element  Formulation 
4.7.1  Strong  form  of  the  problem 

The  governing  field  equations  of  composite  laminates  are  discussed  in  detail  in  Pana- 
handeh,  Masud  and  Ghanimati  [1993].  To  the  general  equilibrium  equations  we  add  the 
boundary  conditions 


Ua 

=  tta 

on 

Fu 

(4.35) 

W3 

=  U3 

on 

r., 

(4.36) 

= 

on 

re 

(4.37) 

In  addition 


=  r“ 

on  F“ 

(4.38) 

=  9 

on  F, 

(4.39) 

=  S° 

on  F, 

(4.40) 

where  PuDFuj  represents  the  boundary’  where  translation  boundary  conditions  are  applied; 
r“  U  r,  represents  to  the  boundary  where  traction  boundary  conditions  are  applied;  Tg 
corresponds  to  the  boundary  with  prescribed  rotation,  and  F,  corresponds  to  the  portion 
of  boundarj-  with  applied  moments. 

As  usual 

n  —  Hq  “  Uq  Aq  (4.41 ) 
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denotes  the  unit  outward  normal  to  the  boundary  d(l  of  domain  Q. 


4.7.2  Weak  form  of  the  problem 

Multiplying  the  strong  form  of  the  problem  with  the  admissible  variations,  integrating 
by  parts  amd  using  the  prescribed  essential  and  the  natural  boundary  conditions  we  arrive 
at  the  weak  form  of  the  problem.  The  spaces  relevamt  to  the  problem  are 


5  =  {(u,U3,ff)j(ti,U3,fi)  e  (u,U3,0)  :  S  -*1l  X  X  71  : 

s.t. «  =  tt  on  Fu,  U3  =  U3  on  Fuj,  on  Fg} 


(4.42) 


V  =  {(u,U3,d)|(«,U3,^)  e  Hq{Q),  {u,U3,9)  :  S  7Z  X  71  X  71  : 
s.t.  u  =  Oon  Fu,  U3  =  Oon  Fu,,  0q  =  Oon  F^} 


(4.43) 


where  5  is  the  space  of  trial  displacements  and  trial  rotations,  and  V  is  the  associated 
space  of  weighting  functions,  respectively.  denotes  the  space  of  square-integrable 

functions  along  with  their  generalized  deri%'atives  defined  over  fl,  and  Hl{Q)  is  the  subset 
of  H^(Q}  wh  ose  members  satisfy  zero  essential  boundary  conditions. 
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4.7.3  Finite  element  matricee  and  vectors 


The  finite  element  stiffness  matrix  and  load  vector  may  be  obtained  directly  from  the 
matrix  form  of  the  variational  equation.  The  finite  element  approximations  for  u,  ti,  B  and 
B  are  denoted  by  and  respectively.  In  a  typical  element,  possessing 

nodes, 

u'*  =  ^  ui  (4.44) 

=  (4.45) 

a=l 

^en 

'■  z=  ^  (4.46) 

0=1 
^en 

'■  =  iV„  Bl  (4.47) 

0  =  1 

where  Na  is  the  shape  function  associated  with  node  ‘a’,  u^,  ti*,  ,  and  are  the  a‘* 
nodal  values  of  ti*,tt*,tf*,and  ^*,  respectively.  It  is  not  necessary  to  assume  that  B^  and 
u*  be  defined  in  terms  of  the  sarnie  shape  functions  and  nodail  patterns. 

4. 7. 3.1  Stress  vectors 

The  resultant  stress  vectors  for  an  element  which  has  the  3D-effects  of  an  elasticity 
element  and  2D  effects  of  a  shell  like  element  are 


(inplane) 


(shear) 
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(4.48) 


(4.49) 
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(bending) 


(through-thickness) 


4. 7.3. 2  Strain  vectors 

The  strain  vectors  corresponding  to  the  stress  vectors  are 


(4.50) 

(4.51) 


(inplane) 


(shear) 


(bending) 


(through-thickness) 


(4.52) 

(4.53) 

(4.54) 

(4.55) 


We  can  combine  equ.  (7. 3. 2.1)  and  (7. 3.2. 4)  to  yield  a  vector  that  incorporates  3D 


effects 


(4.56) 
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4.7.3.3  Matrix  Differential  Operators 


The  strain  vectors  can  be  written  in  terms  of  differential  operators  as  follows 


=  I 


AfdJdX^ 

Ajd/dX^ 

D'^d/dX^ 

Ajd/dX^  +  Ajd/dX^ 


>  u 


(4.57) 


For  the  case  of  flat  geometry  and  linear  analysis,  the  matrix  differential  form  for  the 
modifled  membrane  effect  becomes 


Bm  =  < 


Similarly,  for  the  shear  operator 


djdX^ 

dldX^ 

d/dX^ 

d/dX^  +  d/dX^ 


(4.58) 


~  7q3  —  ~ 


=  [D'^a/dx-  ^1]  {#} 

Therefore  transverse  shear  strain  vector  for  flat  geometries  can  be  speciedized  to 

=  [Bsm  B,fc)  I  “  I 


(4.59) 

(4.60) 


d/dX^  Af 
d/dX^  Aj 


{:} 


(4.61) 


Following  along  the  same  lines,  the  bending  strain  vector  can  be  written  in  terms  of 
matrix  differential  operators 

r 1 

(4.62) 


6o  —  [Bfcni  ^ 

ful 

l^i 

where 

D^d/dX'  1 

Ajd/dX^ 

Bhm  —  * 

Dld/dX^  i  ;B66  =  < 

AldldX'^ 

Dld/dX^  +  Djd/dX'  J 

AfdldX^  +  Ajd/dX^ 

*  t 

>  (4.63) 


Since  for  flat  geometries  D,a  =  0 


d/dx^ 

Btm  =  [0]  ;  Bu  =  d/dX^ 

d/dX^  +  d/dX^ 


(4.64) 


Finally,  a  totsJ  matrix  differential  operator  B  can  be  defined  which  produces  the  total 
strain  vector  i  when  applied  to  the  displacement  field  u  and  the  rotation  field  B: 


where 


-»{:} 


Bm  0 

B,m  Bfb 

Bhm  Bbh 
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4.8  Numerical  Examples 

4.8.1  FVee  edge  boundary-value  problem.  [45,'45]s 

The  first  numerical  simulation  is  a  prismatic  symmetric  laminate  having  traction-free 
edges  at  1'  =  ±a  and  surfaces  Z  —  ±h,  and  loaded  by  strain  applied  only  on  its  ends  at 
X  —  constant.  Each  layer  is  composed  of  unidirectional  fiber- reinforced  materizd  such  that 
the  fiber  direction  is  defined  by  its  angle  6  with  the  X-axis. 

The  elastic  properties  of  various  composite  materials  used  in  this  numericad  simulation 
have  been  talcen  from  Pagano  [1989],  p.  4.  In  this  example  a  laminate  consisting  of  four 
unidirectional  fiberous  composite  layers,  two  with  their  axis  of  elastic  symmetry  (fiber 
direction)  at  +45  and  two  at  —45  to  the  longitudincil  laminate  axis  is  considered.  Figure  4.8 
shows  the  laminate  geometry  and  the  coordinate  system.  1%  strain  in  opposite  directions 
is  applied  at  X  =  —L  and  X  =  L  respectively,  while  it  is  restrained  to  move  in  the  axial, 
lateral  and  thickness  directions  at  X  =  0.  In  order  to  solve  this  problem  a  finite  element 
mesh  comprising  1920  composite  shell  elements  with  2665  nodes  is  generated.  The  physical 
dimensions  for  the  numerical  simulation  are  X  =  60,  T  =  20,  Z  =  2.5,  with  12  elements 
in  X  direction,  40  elements  in  the  Y  direction  and  4  elements  through  the  thickness.  For 
each  finite  element  layer  through  the  thickness,  the  reference  surface  is  assumed  to  be 
associated  with  the  bottom  surface  of  the  composite  shell.  We  will  be  looking  at  a  section 
cut  at  X  =  0  for  the  displacement  and  the  stress  fields. 

Figure  4.9  shows  the  axial  displacement  distribution  at  the  top  free  surface  of  the 
section  cut  at  X  =  0,  auid  the  results  are  compared  with  Moires  et  al.  and  Pagano  et 
al.  [1989].  Figure  4.10  presents  the  interlaminar  shear  stress  at  material  interface  and 
their  corresponding  values  from  J.N.  Reddy  [1993].  Figure  4.11  shows  the  major  stress 
components  <Tii.  <Ti2.  <Ti3,  and  their  corresponding  values  from  the  numerical  simulations 
of  Pipes  et  al.  [1970].  In  order  to  complete  the  discussion,  we  have  2Jso  presented  the 
minor  stress  components  in  Fig.  4.12. 
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4.8.2  FVee  edge  boundary-value  problem  with  a  cylindrical  hole.  [45,-45]s 

The  composite  laminate  considered  in  this  numerical  simulation  has  the  same  physical 
dimensions,  material  properties  and  loading  conditions  as  in  the  previous  case.  However 
the  present  laminate  has  a  \init  diameter  cylindrical  hole  at  (0,0,0)  with  the  axis  of  the 
cylinder  coincident  with  the  Z  axis  (see  Fig.  4.13).  The  ratio  of  the  diameter  of  the  hole  to 
the  width  of  the  laminate  is  0.1.  To  solve  this  problem  a  mesh  containing  720  composite 
elements  per  finite  element  layer  and  with  four  layers  through  the  thickness  was  generated. 
Total  number  of  nodes  in  this  problem  are  3780.  Traction-free  boundary  conditions  are 
applied  on  the  cylindrical  surface  of  the  hole.  The  laminate  is  contrained  to  move  in  the 
axial,  lateral  and  transverse  directions  by  appropriately  constraining  the  nodes  along  the 
symmetry  lines  at  =  0. 

Figure  4.14  presents  the  major  stress  components  trj  j ,  <ri2,  <ri3  at  a  section  cut  at  A  =  0 
in  the  top  layer  with  -1-45  degrees  ply  orientation.  These  results  have  been  compared 
with  the  3D  anisotropic  calculations  done  with  a  mesh  which  is  twice  as  refined  as  the 
composite  mesh.  It  can  be  seen  that  close  to  the  hole  there  is  a  sudden  increase  in  the 
values  of  the  stresses.  In  a  region  away  from  the  free  edge  effects  of  the  traction  free 
boundary  and  the  hole,  the  ratio  of  <rn  and  wij  agrees  closely  with  that  of  the  preceeding 
numerical  simulation.  Figure  4.15  shows  the  minor  stress  components  which  also  show  a 
good  agreement  with  the  results  of  the  anisotropic  3D  simulation. 


Figure  4.8  Laminate  geometry  and  the  coordinate  system. 


Composite  Shell  *  Moires  et.  al  Pagano  et.al 


Figure  4.9  Axial  displacement  distribution  at  the  top  free  surface. 


- Composite  Shell .  J.N.  Reddy 


Figure  4.10  Interlaminar  shear  stress  at  material  interface. 


- SigJI  .  SigJ2  .  SigJ3  Pipes  et.  al 


Figure  4.11  Major  stress  components  in  the  top  -|-45  laminate. 
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Figure  4.14  Major  stress  components  in  the  top  -f45  laminate. 


Sig_^22  - Sig_33  .  Sig_23  Anisotropic 


Figure  4.15  Minor  stress  components  in  the  top  -1-45  l^lminate. 
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4.9.  Conclusions 


In  this  chapter  we  have  presented  a  finite  element  formulation  which  is  suitable  for  the 
analysis  of  multi-director/multi-layered,  shear- deformable  flat  composite  laminates.  The 
geometric  description  employed  for  the  composite  shells  finds  its  roots  in  the  so-called 
degenerated  shell  approach.  A  set  of  kinematic  hypothesis  is  then  introduced  to  accommo¬ 
date  the  effects  of  transversal  warping  of  the  cross  section  due  to  shear  deformation  as  well 
as  fiber  compressibility  resulting  from  transverse  normal  stresses.  The  kinematics  are  in¬ 
dividually  and  independently  represented  for  each  layer,  and  are  coupled  by  the  conditions 
of  displacement  continuity  between  the  layers.  The  rotation  field  is  layer-wise  continuous 
and  is  discontinuous  across  the  layers  in  the  thickness  direction.  IVansversal  warping  of 
the  composite  cross-section  is  described  by  individual  layer  directors  which  simultaneously 
rotate  and  stretch.  This  results  in  discontinuous  strain  fields  across  the  different  mate¬ 
rial  sets,  thereby  creating  the  provision  of  stress  continuity  across  the  material  interfaces. 
Since  the  formulation  resolves  all  strains,  all  stresses  are  computed  through  the  three- 
dimensional  constitutive  equations  and  the  usual  ‘zero  normal  stress’  shell  hypothesis  is 
not  employed.  The  different  material  layers  with  arbitrary  fiber  orientations  may  occupy 
arbitrary  layer  locations  in  the  composite.  In  this  theory,  at  most,  only  first  derivatives  of 
displacement  and  rotation  fields  appear  in  the  variational  equations.  The  practical  conse¬ 
quence  of  this  fact  is  that  only  C°  continuity  of  finite  element  functions  is  required  which 
is  readily  satisfied  by  the  family  of  Lagrange  elements  (Hughes  [1987]).  In  displacement 
formulation  of  plates  and  shells,  special  attention  needs  to  be  given  to  trzinsverse  shear  and 
membrane  terms  to  prevent  the  occurance  of ‘mesh  locking’.  We  have  employed  the  selec¬ 
tive  reduced  integration  technique  on  membrane  and  transverse-sheeir  terms  to  avoid  this 
problem.  Numerical  results  are  presented  to  demonstrate  the  performance  of  the  proposed 
theory. 

The  present  theory  is  suited  for  the  study  of  smart  plate  and  shell  structures  with 
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embedded  sensors  and  actuators.  In  particular  the  use  of  shape  memory  alloys  in  smart 
systems  can  be  modeled  through  the  present  theory  and  the  developments  presented  in 
Chapters  2  and  3. 
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CHAPTER  5 


FUTURE  RESEARCH  AND  POTENTIAL  POST  APPLICATIONS 

It  is  expected  that  wide  spread  application  of  the  intelligent  niaterials  technology  using 
the  present  generation  of  actuators,  sensors,  processors,  and  control  methods  will  occur. 
The  breadth  of  application  of  this  technology  is  expected  to  not  only  include  the  aerospace 
industry  but  to  become  widespread  in  the  home  construction,  automotive,  and  machine 
tool  industries  as  well. 

At  the  present  time,  all  of  the  technologies  needed  for  cost-effective  application  of 
intelligent  structures  have  not  been  sufficiently  developed.  There  are  a  number  of  difficult 
problems  which  remain  to  be  solved.  In  this  section  some  for  the  more  important  of  these 
problems  are  briefly  discussed. 

Actuation  Materials:  In  order  to  truly  achieve  the  desirable  level  of  control  for  many  struc¬ 
tural  applications,  actuation  materials  which  have  3  to  10  times  larger  strain  than  those 
currently  commercially  available  must  be  developed,  Crawley  [1992],  or  materials  should 
be  developed  similar  to  shape  memory  alloys,  but  with  much  higher  bandwidth  than  those 
currently  available.  The  proposed  research  work  will  provide  a  powerful  ansilysis  tool  to 
study  the  optimum  chsu-acteristics  of  these  actuation  materials. 

Development  of  optimization  algorithms  in  order  to  compute  the  most  efficient  sensor 
and  actuators  positions  is  another  objective  for  the  mathematical  modeling  of  adaptive 
structures.  Some  of  these  optimal  placement  strategies  can  be  based  on  finite  element 
analysis  output  such  as  Independent  Modal  Space  Control  which  utilizes  the  eigenvectors 
from  a  finite  element  modal  analysis.  The  proposed  work  has  the  potential  to  be  extended 
to  study  the  optimal  position  of  active  elements  in  intelligent  structures. 

Many  structural  components  in  aerospace  industry  have  curved  geometries.  The  cou- 


pling  of  in-plane  and  out-of-plane  strain  measures  in  these  structures  will  affect  the  in¬ 
teraction  between  sensors/ actuators  and  the  host  materials.  In  particular  for  the  case  of 
shape  memory  alloys  with  large  actuation  strains,  the  results  of  this  research  could  be  used 
to  investigate  the  effect  of  curved  geometries. 

Further  investigation  on  constitutive  modeling  of  shape  memory  alloys  under  cyclic 
thermomechanical  loadings  with  various  preload  conditions  is  needed  for  evaluation  of  the 
integrity  of  smart  structural  systems  in  the  long  term. 

5.1  Potential  Applications 

Incorporating  active  materials  into  the  structure  of  various  aerospace  vehicles  will  al¬ 
low  designers  new  flexibility  with  which  to  increase  the  vehicle  performance.  These  perfor¬ 
mance  enhancements  will  allow  future  aircradt  systems  to  outperform  their  adversaries  and 
better  accomplish  their  missions.  Potential  applications  for  active  structures  span  many 
disciplines  including  aerodynamics,  propulsion,  aeroelasticity,  vibration  and  acoustics  as 
discussed  in  Agnes  &  Silva  [1992]. 

Aerodynamics  applications  for  active  structures  technologies  require  changing  struc¬ 
tural  contours  of  the  mrcraft  to  influence  the  airflow.  Studies  have  shown  that  only  small 
changes  in  airfoil  shape  can  result  in  significant  drag  reduction.  This  application  is  sim¬ 
ilar  to  the  Mission  Adaptive  Wing.  Actuation  requirements  are  for  quasi-static  control 
therefore  shape  memory  alloys  may  prove  best  for  this  application. 

Another  possible  aerodynamic  application  is  the  control  of  unsteady  aerodyn^Lmics  seen 
in  either  high  angle  of  attack  flight  or  in  turbulent  boundary  layers.  This  application  would 
require  local  changes  in  the  airfoil  skin  that  would  interact  with  the  imsteady  airflow,  reduce 
drag,  and  increase  performance. 

Another  set  of  aerodynamic  applications  are  in  propulsion.  Active  flow  field  control 
could  be  used  to  keep  combustion  optimal  in  future  engine  designs.  Active  materials  could 


also  be  utilized  in  thrust  vectoring  nozzles  to  provide  higher  frequency  control  capability 
and  lighter  weight.  Materials  development  will  have  to  proceed  these  applications  due  to 
the  high  temperatures  found  in  aircraft  engines.  The  externals  of  engines,  however,  have 
lower  temperatures,  on  the  order  to  300  —  400°  F,  and  suffer  from  many  vibration  problems, 
particul2irly  in  the  region  of  turbine  rotor  frequencies  (30-100  Hz).  Active  damping  of  the 
external  components  could  provide  life  cycle  cost  savings  for  propulsion  systems. 

Solution  of  the  engine  externals  vibration  environment  would  be  similar  to  other  struc¬ 
tural  dynamics  applications  in  vibration,  acoustics  and  aerr.lasticity.  The  high  vibratory 
environment  of  aircraft  systems  presents  a  range  of  challenges  to  the  designer.  Possible 
applications  would  include  active  isolation  of  stores,  equipment,  sensors,  or  pods;  skin 
panel  fatigue  life  extension;  tail  buffeting  reduction,  active  landing  gear  and  gim-fire  vi¬ 
bration  reduction.  These  and  other  vibration  applications  vary  in  frequency,  bsmdwidth, 
temperature  and  control  authority  requirements.  Solutions  could  involve  active  isolation, 
damping,  control,  or  any  combination  of  the  three. 

Another  class  of  structural  dynamics  applications  is  aeroeleistic.  Coupling  between  the 
aerodynamic  flow  and  structural  dynamics  can  result  in  flutter.  Active  structures  may 
be  able  to  attack  this  problem  from  both  directions:  aerodyneunic  and  structural.  Novel 
control  surfaces  or  active  wing  twisting  may  be  possible.  Direct  active  flutter  control 
through  wing  twisting  will  require  an  increase  in  the  control  authority  currently  available 
from  piezoceramics  or  an  increase  in  the  bandwidth  of  shape  memory  alloys. 

Acoustic  applications  comprise  the  leist  subset  of  structural  dynamic  applications.  Two 
areas  which  have  been  demonstrated  zu'e  cabin  noise  reduction  zind  cavity  resonance  prob¬ 
lems.  In  summation,  a  wide  array  of  aircraft  applications  for  active  structures  technology 
exist;  albeit  at  the  conceptual  stage.  Further  study  in  this  area  is  therefore  required,  and 
this  research  is  a  step  toward  this  direction. 
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